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ABSTPACT 

A  high  resolution  acoustic  navigation  system  for  ocean  use 
is  being  developed  at  the  Woods  Hole  Oceanographic  Institution. 
The  system  can  yield  navigation  fixes  with  respect  to  a  bottom 
moored  reference  net  with  accuracies  (on  a  fix  to  fix  basis)  of 
a  few  centimeters.   In  order  to  use  the  system  to  best  advan- 
tage a  survey  is  required  to  determine  precisely  the  relative 
positions  of  the  net  elements.   Each  element  combines  a  pulse 
transponder  with  a  continuous  wave  (CW)  beacon.   Accumulated 
phase  (Doppler  shift  of  the  CW  beacon)  between  survey  points  is 
measured  as  well  as  acoustic  travel  times  between  survey  points 
and  transponders.   Non-linear  regression  techniques  are  employed 
to  develop  a  maximum  likelihood  estimator  for  net  element 
positions  based  on  these  phase  and  travel  time  measurements. 
An  approximate  error  covariance  matrix  is  generated  and  an 
optimum  choice  of  survey  points  is  indicated.   The  combined 
system,  using  these  selected  locations  for  performing  the  survey, 
can  yield  reference  mooring  coordinates  with  error  of  ±1  meter. 
Improved  precision  appears  to  be  limited  by  inaccuracies  in  the 
pulse  and  Doppler  measuring  system. 
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I  INTRODUCTION 

There  are  many  methods  currently  used  to  position  ships, 
submersibles,  buoys  or  submerged  instruments  in  the  deep  ocean. 
Broadly  speaking  they  can  be  divided  into  two  types.   The  first 
method  positions  via  electromagnetic  transmission  from  a  shore 
station  or  satellite  whose  position  or  orbit  is  precisely  known. 
Such  systems  operate  at  long  ranges  but  are  limited  in  accuracy 
to  about  100-200  meters  at  best.   Furthermore,  these  systems 
require  above  surface  antennas  and  are  therefore  of  little  use 
in  positioning  entirely  submerged  vehicles.   The  second  method 
positions  via  acoustic  transmission  from  a  set  of  transponders 
or  beacons,  usually  moored  to  the  ocean  bottom,  whose  relative 
positions  are  precisely  known.   When  operated  in  the  pulsed  or 
transponder  mode,  position  inaccuracies  of  an  acoustic  system 
can  be  as  small  as  10-20  meters  in  5  km  water  depths.   When 
operated  in  a  continuous,  or  Doppler  mode,  the  accuracy  can 
approach  3-4  centimeters.   The  range  of  acoustic  systems  is 
usually  severely  limited  by  a  number  of  factors,  the  most  im- 
portant of  which  are  acoustic  refraction  and  attenuation  of 
sound  in  sea  water  due  to  spreading  and  absorption  in  the  ocean. 
Current  acoustic  systems  have  ranges  of  about  10  km. 

A  navigation  system  utilizing  acoustic  beacons  (Fig.  1) 
which  operate  in  both  a  pulse  and  continuous  wave  mode  is  being 
developed.   This  system  capitalizes  on  the  attributes  of  both 
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TOWED  HYDROPHONE 


Figure  1.   Acoustic  Navigation  Beacon  Net 
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the  pulse  and  Doppler  modes  and  will  be  capable  of  positioning 
a  platform  with  respect  to  an  array  of  reference  beacons  to  an 
accuracy  of  1-2  meters  and  of  repositioning  a  platform  within 
10  centimeters  of  a  previous  position.   To  make  full  use  of  the 
system  capabilities  a  survey  is  required  to  accurately  deter- 
mine the  relative  positions  of  the  reference  beacons.   The  pur- 
pose of  this  investigation  is  to  develop  an  accurate  survey 
technique.   Algorithms  for  estimating  the  beacon  positions  and 
the  errors  in  those  positions  are  derived.   Performance  is 
predicted  using  computer  models.   The  estimation  method  used  is 
a  straightforward  extension  of  the  pulse  system  technique  des- 
cribed by  Hunt  et  a_l.  (1974)  . 

1)  Description  of  System 

Pulse  (transponder)  navigation  systems  measure  the  travel 
time  of  an  acoustic  pulse  to  estimate  the  slant  range  to  a 
transponder.   Shipboard  processing  equipment  interrogates  the 
bottom  moored  transponders  and  measures  the  time  of  arrival  of 
the  reply  pulse  from  each  transponder.   Sound  velocity  correc- 
tions and  corrections  for  fixed  system  delay  times  are  applied, 
and  travel  times  are  converted  into  estimates  of  slant  range  to 
the  transponders.   The  accuracy  of  the  measurement  is  depen- 
dent upon  the  pulse  bandwidth  and  the  signal-to-noise  ratio  in 
the  received  pulse.   Systems  of  this  type  are  widely  in  use  in 
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oceanographic  research,  submerged  rescue  operations  and  commer- 
cial enterprises  in  the  ocean  (Boegeman  et.  a_l.  ,  1972;  Cestone 
and  St.  George,  1972;  Van  Calcar,  1969). 

7A  continuous  wave  (CW)  beacon  system  has  been  developed 
at  the  Woods  Hole  Oceanographic  Institution  to  accurately 
track  the  position  of  ship-mounted  and  submerged  hydrophones 
(Porter  et  al_.  ,  1973).   The  system  measures  the  frequency 
change  of  a  continuous  tone  due  to  hydrophone  motion.   Hydro- 
phone velocity  and  displacement  are  derived  from  the  beacon 
signal  Doppler  shift  which  is  given  by 

Af  =  fB  -  fD  =  -fB(v/c), 

where  fD  is  the  Doppler  shifted  beacon  frequency  due  to  a 

hydrophone  velocity  v,  f   is  the  beacon  frequency  and  c  is  the 

b 

speed  of  sound  in  water.   Doppler  shift  is  proportional  to 
beacon  frequency.   Hence  the  accuracy  of  the  measurement  is 
dependent  on  the  beacon  frequency  and  the  signal-to-noise 
ratio  in  the  narrow  band  that  encompasses  the  maximum  Doppler 
shift.   Because  the  beacon  frequency  is  much  greater  than  the 
pulse  bandwidth  of  transponder  systems  the  Doppler  system  can 
measure  range  changes  much  more  accurately  than  the  pulse 
system  assuming  that  both  systems  have  the  same  signal-to- 
noise  ratio.   In  the  present,  implementation  of  the  system 
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accumulated  phase  change  0 (t)  =  2  ir(Af)t  is  measured  by  digital 
techniques  and  hydrophone  displacements  of  a  quarter  wavelength 
can  be  resolved.   A  wavelength  for  13  kHz  (a  typical  beacon 
frequency)  is  /\   =  11.5  centimeters,  so  that  resolution  is 
about  X   /4  ~  3  cm.   Phase  changes  due  to  inhomogeneities  and 
fluctuations  in  the  water  column  are  negligible  over  the  range 
of  the  system.   We  discuss  this  in  Section  III.  2. 

2)  Previous  Survey  Techniques 

Several  investigations  have  been  conducted  to  develop  an 
optimum  technique  for  determining  the  coordinates  of  ocean 
bottom  acoustic  transponders.   No  satisfactory  survey  has  been 
developed  previously  for  an  array  of  bottom-moored  CW  beacons. 

Surveys  of  acoustic  transponders  generally  fall  into  two 
categories:  baseline  crossing  and  iterative  techniques.   The 
baseline  crossing  technique,  described  by  Haehnle  (1967)  and 
Hart  (1967),  is  often  referred  to  as  the  conventional  or  class- 
ical transponder  survey  method.   The  method  requires  that  the 
depth  of  the  transponders  be  previously  determined  from  inde- 
pendent measurements.   The  baseline  length,  i.e.,  the  horizon- 
tal distance  between  two  transponders/  is  determined  by  steaming 
across  a  baseline  and  measuring  the  slant  range  to  the  two 
transponders  during  the  traverse.   This  requires  accurate  ship 
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positioning  and  large  amounts  of  ship  time  to  yield  accurate 
results.   More  recent  iterative  techniques  are  based  on  or  are 
similar  to  that  described  by  Vanderkulk  (1961).   Vanderkulk 
showed  that  for  three  bottom  transponders  and  six  co-planar 
survey  positions  whose  depths  are  known,  the  positions  and 
depths  of  the  transponders  can  be  found  by  solving  six  linear 
equations.   A  unique  solution  exists  when  the  survey  positions 
do  not  lie  on  a  conic  section.   When  more  than  six  survey  posi- 
tions are  used  the  additional  equations  result  in  an  over- 
determined  situation  for  which  minimum  least  squared  error 
techniques  have  been  successfully  applied,  (Lowenstein,  1965; 
McKeown,  1974;  Mourad  et  a_l .  ,  1972;  Heckman  and  Abbott,  1973; 
Hunt  et  a_l.  ,  1974)  . 

3)  Pulse-Doppler  Survey  Technique 

The  pulse-Doppler  survey  consists  of  making  measurements 
of  travel  times  of  acoustic  pulses  and  accumulated  phase 
changes  of  CW  signals  at  various  survey  points.   These  are  con- 
verted to  slant  ranges  and  slant  range  differences,  respect- 
ively.  The  pulse-Doppler  survey  can  yield  more  accurate  results 
than  the  pulse  survey  alone  due  to  the  higher  resolution  of  the 
Doppler  system  measurement. 

The  slant  range  (SR)  between  a  transponder  and  a  survey 
point  is  equal  to  the  travel  time  (TT)  multiplied  by  an 

-12- 


effective  sound  speed,  Ve : 

SR  =  Ve  •  TT 
The  effective  sound  speed  is  a  function  of  the  sound  velocity 
profile,  the  depths  of  the  transponder, and  the  survey  point 
and  varies  with  the  amount  of  acoustic  refraction  (Vass,  1966) . 
For  m  survey  points  there  are  m  measurements  of  slant  range 
for  each  transponder. 

The  change  in  slant  range,  DSR,  between  a  beacon  and  two 
survey  points  is  proportional  to  the  change  in  accumulated 
phase,  0,(t2)  -  ^(t,  )  and  the  wavelength  of  the  signal,  \    -   c/f 

DSR  =  [  0(t2)    -   0(tL)]   X  /2  tt 
The  received  signal, 

s(t)  =  A(t)  exp  {  -j2  it  f  .t  }  , 
where  f~.   =    fB  +fBv/c  is  the  Doppler  shifted  frequency,  can 
be  written  as 

s(t)  =  A(t)  exp  {  -j  [2  tt  ft  +  0(t)   ]  } 
where  the  phase  variation  due  to  a  finite  receiver  velocity  is 
0{t)    =  2  tt  ffit  v/c  (Porter  et  a^L.  ,  1973).   In  general  v  is  a 
function  of  time,  v  =  v(t).   However,  we  assume  that  the 
sampling  interval  is  short  enough  to  consider  v  constant.   in 
practice,  an  interval  of  .1  sec  is  used  so  the  assumption  is 
valid.   Measurement  of  accumulated  phase  change  provides 
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velocity  and  displacement  information  along  the  beacon-survey 
point  vector,  i.e.,  the  displacement  is  the  change  in  slant 
range.   For  m  survey  points  there  are  m  -  1  measurements  of 
slant  range  differences  for  each  beacon. 

The  measurements  of  travel  times  and  phase  changes  may 
be  made  at  the  same  or  different  survey  points.   Only  the  case 
where  the  measurements  are  made  at  the  same  points  will  be 
discussed  in  detail  since  this  technique  is  more  efficient  in 
practice.   The  method  of  solution  for  the  transponder  coordin- 
ates is  the  same  in  either  case  but  the  associated  errors  will 
be  different. 

Non-linear  regression  methods  as  outlined  by  Draper  and 
Smith  (1966)  are  used  to  analyze  this  survey  problem.   Under 
the  assumptions  given  in  the  following  section  the  maximum 
likelihood  estimator  for  the  beacon  positions  is  found.   An 
approximate  covariance  matrix  for  the  estimates  and  a  pro- 
cedure for  optimizing  the  choice  of  survey  points  is  developed 
The  sources  of  error  and  their  effect  on  the  estimates  of  the 
beacon  coordinates  are  discussed. 
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II  PARAMETER  ESTIMATION 

The  problem  of  parameter  estimation  can  be  briefly  des- 
cribed as  follows:  given  a  system  from  which  information  can 
be  received  and  whose  state  is  characterized  by  a  set  of  p 
parameters,  make  some  estimate  of  the  state  of  the  system  from 
the  information  received.   The  information,  which  will  gener- 
ally be  perturbed  by  noise,  may  be  measurements,  a  message,  or 
some  other  observations  on  the  system. 

The  following  definitions  are  used: 

0/  =  [9^,  92/.../©p]   =  true  parameter  values 

8/  =  [9,,  9^ , . . . , 9  ]   =  estimated  parameter  values 

F'  =  F_'  (9)  =   [f  •,  ,  f2  ,  .  .  .  f  ]  =  noiseless  message 

Y'  =  [y-i /Y2  '  •  •  •  »  yn  ]   =  messa9e  (observations  perturbed 

by  noise) 
where  9_,  9_,  F_,  and  Y  denote  column  vectors  and  (  '  )  denotes 
transpose  (Manasse,  1960) . 

When  a  noiseless  message  F(j|)is  received,  the  information 
is  sufficient  to  determine  9_  exactly  if  the  number  of  measure- 
ments n  is  not  less  than  the  number  of  parameters  p.   When 
n  =  p  the  problem  is  called  the  MINIMUM  DATA  case.   When  n  >  p 
it  is  called  the  REDUNDANT  DATA  case. 
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1)  The  Maximum  Likelihood  Estimate 

When  the  observations  Y  are  received,  some  criteria  must 
be  used  to  determine  an  optimum  estimate  Q_   since  Y  differs 
from  F_  due  to  the  presence  of  noise.   One  often  used  decision 
rule  is  based  on  the  maximum  likelihood  principle.   This 

A 

principle  prescribes  that  the  observer  choose  the  9_  which 
renders  the  observations  Y  most  likely;  that  is,  it  maximizes 
the  conditional  probability  density  P  ( Y  |  9_)  .   When  expressed 
as  a  function  of  Y  and  9_  this  is  called  the  likelihood  function 
and  is  sometimes  written  L(Y;9).   When  the  a_  priori  probability 
density  P  ( 9_)  can  be  considered  to  be  a  constant  over  the 
region  of  interest  then  the  maximum  likelihood  estimator  also 
maximizes  P(9_|y).   This  can  be  shown  by  the  use  of  Bayes 
Theorem  to  express  the  a_  posteriori  probability  density 
P  (9  |  Y)  in  terms  of  P(Y  [  9_)  and  P(9)  : 


P(Y|9)  P(9) 

P(9|Y)  =  =  C1P(Y|9)  (1) 

P(Y) 


where        P  ( Y)  =  /p(y|£)  p  (Ji)  d^. 

(Note  that  P  ( Y)  is  not  a  function  of  9^  and  can  be  treated  as  a 
constant  of  proportionality) . 

The  estimate  9_  which  maximizes  the  likelihood  function 
L(Y;9_)  must  satisfy  the  equations 
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Small  variance  is  a  desired  property.   An  asymptotic  measure 
of  this  property  is  efficiency.   An  estimator  is  said  to  be 
efficient  if  the  asymptotic  variance  of  the  estimate  is  no 
larger  than  the  asymptotic  variance  obtained  using  any  other 
estimator.   The  maximum  likelihood  estimator,  when  it  is  con- 
tinuous and  the  first,  and  second  derivatives  exist  and  are 
absolutely  integrable,  is  efficient.   It  can  be  shown  that 
maximum  likelihood  estimates  are  asymptotically  normally  dis- 
tributed  with  mean  9  and  variance  1/R  (9)  where 


R2(S)=  E 


/  ^log  L 


(Kendall  and  Stuart,  1967), 


89 

In  the  case  of  additive  Gaussian  noise  the  form  of  the 
likelihood  function  can  be  determined  and  expressions  for  the 
variance  of  the  estimates  can  be  developed. 

3)  The  Maximum  Likelihood  Estimate  in  the  Case  of  Additive 
Gaussian  Noise 

The  message  is  assumed  to  be  perturbed  by  additive 
Gaussian  noise,  e,    so  that 

Y  =  F(9)  +  e  (2) 

where  e_  is  characterized  by  an  n  x  n  moment  matrix  M; 

M=E(   [Y-F][Y-F]')  =E  (ee  '  ) 
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The  multidimensional  Gaussian  distribution  of  the  error  can 
be  expressed  as 


P(e)  =  C2  exp  f  -1/2  [e'M^e]  ]• 


or 


P(Y  -  F)  =  C2  expl-l/2[Y  -  F  'M_1  'y  -  F_!  1 
where  C2  is  a  positive  constant  (see,  for  example,  p.  151 
Davenport  and  Root,  19  58) . 

From  equation  (2)  it  can  be  seen  that  the  probability 
density  of  Y  given  F_  is  equal  to  the  probability  density  of 
the  error;  that  is 

P(Y|F)   =  P(e)  =  P(Y  -  F) 
and  since  F  is  a  function  of  9_  this  density  is  also  the  proba- 
bility density  of  Y  given  Q.       The  likelihood  function,  there- 
fore, is  equal  to  the  probability  density  function  of  the 
error,  namely  a  multidimensional  normal  distribution,  that 
is, 

L(Y;9)  =  P(Y|0)  =  P(Y  -  P (9) ) . 
Thus 

L(Y;9)  =  C2  exp  {  -1/2  [y  -  FT  M_1  [y  -  F  J  |  (3) 
The  likelihood  function  is  a  maximum  when  the  magnitude  of 
the  exponent,  |~Y  -  F_T  M~^  l"Y  -  Fl  ,  is  a  minimum. 

When  the  noise  e_  has  zero  mean  and  the  error  in  each 
measurement  is  statistically  independent  of  the  errors  in  the 
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other  measurements  the  moment  matrix  M  is  diagonal 

.2 
M  =  E(ee' ) 


°i  a-*     0 


0  <f2 

n 


and  M 


-1 


0 


•i/cr 

n 


.  2  -t-v> 

where  OT   denotes  the  variance  of  the  error  in  the  k 
k 


measure- 


ment.  The  moment  matrix,  M,  can  be  written  in  terms  of  weight- 
ing factors  w-  and  the  standard  error  u  .   The  matrix  of 


weighting  factors  is  denoted  by  V 


-1 


M  =  vcT 


■/ 


and 


w- 
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w- 


0 


w2   0 
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(4) 


2     2 

where  w.  =  (f   /  (J      . 

1  i 


The  maximum  likelihood  estimate  in  this  case  is  the  weighted 
least  squares  estimate  since  the  magnitude  of  the  exponent 
in  equation  (3)  is  proportional  to  the  weighted  sum  of  squared 
errors: 
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n 


i  y 

1/2  [X  -  F(9)]'  M~1[y  -  F(9)J  =  — g-  /^  '  (Yi  -  fi(9))2  W± 
When  the  variances  of  the  errors  are  equal  then  the  weighting 
factors,  w^,  are  equal  to  one,  and  the  maximum  likelihood 
estimate  is  the  familiar  least  squares  estimate  which  minimizes 
the  sum  of  squared  errors  SS(9_): 

n  2 

SS(9)  =   £    (y±  -  fi(9))  ■ 

i=l 

Assume  that  a  rather  accurate  estimate  of  9^  has  been 
obtained  and  that  F_(9_)  can  be  approximated  by  a  first  order 
Taylor  series 

F(9)  =  F(9)  +  X  A9  , 

9  9    9 


where  X 


is  an  n  x  p  matrix  of  derivatives  with  ij 


element   <3f  •       and  the  errors  Z_,  9  =  9_  -  9_  are  assum-d  small. 

When  this  expression  for  F_(9_)  is  substituted  in  equation  (3), 
and  when  the  a_  priori  probability  density  P  (9_)  can  be  consider- 
ed to  be  a  constant  over  the  region  of  interest,  the  probability 
density  of  9_  given  Y  can  be  shown  to  be  a  multidimensional 
normal  distribution  with  moment  matrix    : 


P(9|Y)  =  C3  exp^2  A9'  n'Ae} 


(5) 
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n"1      -l    -2    -i 

where   _>_    =  X'M  X  =  g-      (X'V   X)  .  (6) 

The  moment  matrix  J is  known  as  the  variance-covariance 

matrix  since  it  consists  of  the  variances  (diagonal  terms)  and 
covariances  (off-diagonal  terms)  of  the  estimates.   it  is  also 
known  as  Fisher's  information  matrix  (Van  Trees,  1968).   Specif- 
ically, under  the  conditions  of  uncorrelated,  zero  mean  Gaussian 
noise,  the  variance  of  the  estimates  can  be  expressed  as 

2^2 


wi      'n 


-2 


I 

i=l 


w 


3  9./     1 
l  / 


-1 


This  estimate  is  approximately  equal  to  the  asymptotic 

2 

variance  1/R  (9)  is  close  to  9. 


4)  Application  of  Parameter  Estimation  Techniques  to  the  Survey 
Problem 

The  pulse-Doppler  navigation  system  makes  measurements 
which  are  characterized  by  a  set  of  unknown  parameters.   The 
unknown  parameters  are  the  geometrical  coordinates  of  the 
beacons  and  the  ship  positions  at  which  the  measurements  are 
taken.   The  measurements  consist  of  slant  ranges  and  slant 
range  differences. 

An  orthogonal  XYZ  coordinate  system  will  be  used  in  which 
the  beacon  coordinates  are  (0,  0,  Z..  )  ,  .  (X~  ,  0,  Z2  )  ,  and  (X-^Y-^Z-) 
as  shown  in  Figure  2.   (Only  the  3  beacon  survey  problem  will 
be  considered  here;  the  method  can  easily  be  extended  to  more 
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than  3  beacons) .   The  depth  of  the  ship  hydrophone  is  assumed 
to  be  a  known  constant,  ZS.   Simultaneous  measurements  of  slant 
ranges  and  slant  range  differences  are  taken  at  m  survey  points 
(XSi/YSi,ZS) . 

The  unknown  parameters,  that  is,  the  beacon  and  ship  posi- 
tion coordinates,  can  be  expressed  as  a  column  vector  0_,  with 
transpose 

6   =  [0,,0p,...,0  J 

where  p  =  2m  +  6. 
The  "noiseless  measurements"  are  the  exact  geometrical  slant 
ranges  and  slant  range  differences.   These  can  be  expressed  as 
SRij  =  SRij(t)  =   [(XSi-Xj)2  +  (YSi-Yj)2  +  (ZS-Zj)2  ]1/2 
where  X-L  =  Y1  =  Y2=0     i  =  l,2,...,m 

j  =  1,2,3 
Similarly,  the  slant  range  differences  between  two  survey  points 

i  and  ii  and  the  j    beacon  can  be  expressed  as 

2  2  ?    1  /2 

DSRijte)  =   [(XSi;L-Xj)   +  (YSii-Yj)   +  (ZS-ZjT  ] 

2  2  2    1/2 

-  [  (XS^Xj.)   +  (YS^YjT  +  (ZS-Zj)   ] 

These  functions  can  be  expressed  as  a  column  vector  F(9_) 
with  transpose, 

F'  (0)  =   [  SR11,SR12'SR13'SR21'  '  *  * 'DSRH/DSR12,DSR13/  *  *  *  "*  * 
For  m  survey  points  F(9_)  has  dimensions  3m  +  3  (m-1)  =  3(2m-l). 

If  F(0)  could  be  measured  without  errors  then  some  minimum 
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number  of  measurements  would  be  sufficient" to  determine  0_ 
exactly.   In  this  case  the  slant  range  differences  are  linear 
combinations  of  the  slant  ranges  and  the  minimum  data  case,  as 
shown  by  Vanderkulk  (1961),  requires  18  slant  range  measurements: 

SR±j (0) ,   i  =  1,2,. ..,6,   j  =  1,2,3 
or  equivalently,  3  slant  range  measurements  and  15  slant  range 
differences 

SR-lj  (0)  ,   j  =  1,2,3 

DSR..(0),   i=l,2,...,5,   j  =  1,2,3   . 
If  only  slant  range  differences  can  be  measured  then  24  measure- 
ments, 

DSR..(0),   i  =  1,2,.. .,8,    j  =  1,2,3 
are  sufficient  to  determine  0_  exactly  provided  one  has  reasonable 
initial  estimates  of  the  parameters.   This  condition  is  a  result 
of  the  linear  approximation  used  to  solve  the  non-linear  equations 

The  observations  Y  differ  from  F(0_)  due  to  the  presence 
of  noise.   Slant  range  observations  from  the  transponder  system 
can  be  expressed  as 

sij  =  SRij  (i)  +  e  ij 

where  errors,   C- •»  are  assumed  to  be  independent  identically 
distributed  normal  random  variables  with  variance  (X  .   obser- 
vations of  slant  range  differences  from  the  Doppler  system  can 
be  expressed  as 

DS. .  =  DSR. • (0)  +  §.  . 
ID       1J         ID 
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where  the  errors,  O  .  .,    are  assumed  to  be  independent  identi- 
cally  distributed  normal  random  variables  with  variance  (fc    . 
The  errors  for  the  complete  system  can  be  expressed  as 
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[y  -  F(9)]  . 


If  the  errors  £  and  u   are  independent  of  each  othe^  than  e_ 
is  characterized  by  the  moment  matrix 


M  =  VcT  = 


tfiV    o 

0   '  •  -On2 


rt  -  qr2 


i  =  1,2,..., 3m 


tf3m+i  =  0^2 


i  =  1,2, ... ,3 (m-1) 

5)  Gauss-Newton  Method  of  Non-linear  Least  Squares  Estimation 
We  have  seen  that  the  maximum  likelihood  estimator  under 
conditions  of  additive,  uncorrelated,  Gaussian  noise  with  zero 
mean  is  equivalent  to  the  weighted  least  squares  estimator  and 
that  this  estimator  is  consistent  and  efficient.   The  weighted 

A. 

least  squares  estimate  is  the  0^  which  minimizes  the  residual 


sum  of  squares 

SS(9)  =  e'V"1^  =  [y  -  F(9) 


V 


-1 


Y  -  F(9) 


(7) 


In  general  the  estimate  9_  is  the  solution  of  the  p  normal 


equations 


9  SS (0 ) 
d  9i 


=  o   i  =  1,2, . . . ,p 


assuming  that  the  solution  is  interior  to  the  parameter  space. 
When  these  equations  are  linear  they  can  be  solved  directly. 
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When  the  normal  equations  are  non-linear  as  they  are  in  this  case 
an  iterative  technique  must  be  used.   A  number  of  techniques 
exist;  for  example,  quasi-Newton,  conjugate  gradient,  etc.   A 
brief  review  of  numerical  techniques  for  fitting  non-linear 
models  is  given  by  Chambers  (1973). 

a 

The  method  used  to  find  the  estimate  0_  is  based  on  a  modi- 
fied Gauss-Newton  procedure.   The  Gauss-Newton  procedure  consists 
of  linearly  approximating  the  function  F_(9J  and  applying  standard 
linear  least  squares  techniques. 

One  linear  approximation  of  F(9_)  about  an  initial  estimate 
9n  is  a  simple  Taylor  series 

F(9)  =  F(9q)  +  X  A9 

is  an  n  x  p  matrix  of  derivatives  with 


A 


3  F 

A 

• 

wnere  a  — 

-         d  e 

1  s 

elements   x—    = 

3fi 

36  j 

A 

and  A0  =  0  -  0   is  a  column  vector 

-   -o 


of  length  p.    Choosing  the  new  variable  Z_  =  Y  -  F(9_q)  we  have 
the  linear  model 

z  =  x  Ae  +  e. 
The  linear  approximation  to  the  residual  sum  of  squares  is 
SS(9)  =  e'V"1^  =  fz-X  A  9  1'V-1  f  Z -X  Ael=  SSCBq)  , 
The  Gauss-Newton  step  Aq  is  obtained  by  linear  regress- 
ion techniques  on  this  linear  residual  sum  of  squares 
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(Draper  and  Smith,  1966): 

Ae  =  (x,m""1x)~1  x'y^z 

A  new  estimate  0,  =  9_q  +A9  can  be  formed  and  the  process 
repeated,  until  a  given  convergence  criteria  is  reached. 

6)  Modifications  of  the  Basic  Gauss-Newton  Procedure 

Modifications  of  the  basic  Gauss-Newton  procedure  can  be 
made  to  insure  that  the  step  A  9,  will  reduce  the  residual  sum 
of  squares  and  avoid  singularity  problems  in  the  X'M  X 
matrix.  A    complete  discussion  of  the  modifications  used/ 
based  on  stepwise  regression  techniques,  is  contained  in 
Jennrich  and  Sampson  (1967).   A  summary  of  that  discussion  is 
presented  here. 

One  important  modification  of  the  basic  Gauss-Newton 
procedure  is  the  use  of  partial  steps,  7?  A  9,  in  place  of  the 
full  step  A  9  when  the  full  step  results  in  an  increase  in  the 

A 

residual  sum  of  squares  SS(9_).   This  frequently  occurs  when 
the  linear  approximations  upon  which  the  method  is  based  fail. 
The  proportion  W  to  be  used  is  obtained  by  trying  the  arbitrary 
sequence  of  values  1,1/2,  1/4,  1/8,...  until  the  residual  sum 
of  squares  is  reduced.   It  can  be  shown  that  a  sufficiently 
small  step  in  the  direction  of  A  9  will  always  reduce  the 

A. 

residual  sum  of  squares  unless  9_  is  already  a  minimizing  value. 
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Another  modification  of  the  basic  procedure  is  dictated 
by  the  possibility  that  the  X'X~  2£  matrix  is,  or  is  nearly, 
singular.   Step-wise  regression  techniques  are  used  where  the 
parameter  selected  at  a  given  step  is  the  one  which  makes  the 
greatest  reduction  in  the  residual  sum  of  squares.   if  there 
is  a  singularity  problem  then  only  a  subset  of  the  parameters 
is  used.   This  also  provides  a  convenient  way  to  handle  bound- 
ary restrictions  on  the  parameter  values.   A  parameter  is 
modified  only  if  the  new  parameter  value  is  within  or  on  the 
parameter  boundary. 

7)  Iterative  Technique 

In  summary  the  procedure  used  consists  of  the  following 
steps: 

1)  The  value  of  the  function  F_(9_)  and  its  derivatives  X 
are  calculated  using  an  initial  estimate  Q_q. 

2)  The  X'V   X  matrix  is  formed  and  the  residual  sum  of 
squares  SS(9^)  is  computed. 

3)  The  X'V  X  matrix  is  inverted  in  a  step-wise  manner 
and  the  step  size  A   9  =  (X'V  X)"  2£'X_1Z  is  formed. 

4)  A  new  estimate  9i  =  9n  +  A  9   is  formed  and  the  new 
residual  sum  of  squares  SS  ( 9,  )  computed.   If  SS(9_^)  is  larger 
than  SS(90)  then  the  step  size  is  halved  until  SSfe,)  -    SS(9Q) 


o 


10 


r  until  7J  =  (1/2)  iU. 
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5)  This  completes  an  iteration.   The  process  is  repeated 
until  the  solution  converges.   Convergence  is  determined  by 
the  relative  change  in  the  residual  sum  of  squares,  i.e.  the 
solution  has  converged  when  the  relative  change  is  less  than  a 
specified  criterion,  C, 

I  ss(®n+l>  ~  SS^n)| 

' =  C  '  <  C  . 

ss(en+1) 

6)  In  practice  we  require  that  4  successive  iterations 
result  in  C'<  C  in  order  to  accept  the  hypothesis  of  conver- 
gence.  In  some  cases  the  algorithm  may  fail  to  converge  to  a 
minimizing  solution  (e.g.  when  it  converges  to  a  local,  rather 
than  a  global,  minimum) .   A  plot  of  the  parameter  estimates 
and  examination  of  the  residuals  will  usually  provide  some 
indication  of  this. 

8)  Error  Estimates 

We  have  seen  that  maximum  likelihood  estimates  are 
asymptotically  normally  distributed  with  known  variances  and 
that  these  variances  are  the  diagonal . elements  of  the  variance- 
covariance  matrix  |_  .   For  a  process  perturbed  by  uncorrelated 
Gaussian  noise  with  zero  mean  this  matrix  is  given  by  the 
expression 


.  =  (X'V-1^)  1(/2  where  V  1 


W-i 

w2 

0 

o-' 

•w 
n 
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The  (X'V  X)  matrix  is  called  the  correlation  matrix  and  its 

inverse  will  be  denoted  by  A. 

The  standard  error  in  estimating  the  i    parameter,  0 

y . 

1 

is  the  square  root  of  the  i    diagonal  element  of  the  variance- 
covariance  matrix    _: 

(8) 


ei        ! 


where  a.  is  the  i    diagonal  element  of  the  A  matrix 


n 

£ 

k=1 


3f) 

3e 


W, 


(9) 


J 


The  terms  on  the  right  hand  side  of  equation  (8)  are  indepen- 
dent.  The  standard  error,  u,    depends  upon  the  accuracy  with 
which  the  slant  ranges  and  slant  range  differences  can  be 
measured.   An  estimate  of  the  standard  error  can  be  obtained 
from  the  residual  sum  of  squares 


o-2  = 


n-p 


SS(9) 


where  SS(9)  is  defined  in  equation  (7). 

The  other  term,   V  a.,  is  sometimes  called  the  error 

l 

magnification  factor  and  is  a  function  of  the  survey  geometry, 
The  magnitude  of  the  error  magnification  can  be  obtained  by 
evaluating  equation  (9) . 
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The  approximate  errors  in  estimating  the  beacon  coordin- 
ates 

(7q   =  <f^~*\       i   =  1/2,...  ,6 
^i 

can  be  reduced  in  two  ways:  1)  improvements  in  the  measurement 
accuracy  and  2)  a  suitable  selection  of  the  survey  coordinates 
In  the  next  section,  we  will  consider  in  detail  how  this  may 
be  accomplished. 
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Ill  MINIMIZATION  OF  ERRORS 

1)  Measurement  errors  (pulse  system) 

Slant  ranges  and  slant  range  differences  are  measured  by 
two  different  systems  for  which  the  magnitude  and  sources  of 
error  will  in  general  be  different.   A  discussion  of  these 
errors  requires  a  complete  description  of  the  system  operations. 

The  pulse  system  is  designed  to  measure  round  trip  travel 
time  of  an  acoustic  pulse  between  a  ship-mounted  or  ship- 
suspended  transducer  and  near-bottom  transponders.   When  an 
acoustic  pulse  is  transmitted  from  the  transducer,  digital 
counters  in  a  shipboard  receiver  begin  measuring  elapsed  time. 
The  transponders  receive  and  recognize  the  transmitted  pulse 
and  generate  a  "reply"  pulse  at  a  specific  frequency.   When 
a  return  pulse  is  detected  by  the  shipboard  receiver  at  a 
given  frequency  the  corresponding  counter  is  stopped  and  the 
elapsed  time  is  displayed  and  transferred  to  a  digital  com- 
puter for  processing. 

There  are  fixed  time  delays  associated  with  this  process; 
e.g.  signal  recognition  time,  transponder  turn  around  time 
and  signal  processing  time.   Accurate  calibration  of  the  system 
can  reduce  the  uncertainties  in  each  of  these  delay  times  to 
approximately  ±.5  milliseconds  (msec).   For  most  applications 
these  timing  errors  can  collectively  be  considered  to  be  normally 
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distributed  with  zero  mean  and  standard  error  of  1  msec.   Two 
of  these  timing  errors,  signal  recognition  time  at  the  trans- 
ponder and  at  the  ship  are  functions  of  the  signal-to-noise 
ratio  and  the  pulse  bandwidth.   For  example,  the  errors  in 
recognition  time  increase  to  approximately  .8  to  1  msec  for  a 
signal-to-noise  ratio  at  the  pulse  receiver  of  27  dB.   This 
signal-to-noise  is  typical  of  a  range  of  10  km  in  sea  state  3 
for  the  system  presently  in  use  (source  level:  189  dB,  re: 

I  p Pa  pulse  length:  10  msec).   For  ranges  greater  than  about 

II  km  (for  the  present  system  configuration)  signal  recognition 
is  erratic  and  the  errors  can  be  assumed  infinite.   For  fixed 
bandwidth,  timing  errors  can  be  reduced  by  increasing  signal 
power  or  they  may  be  acounted  for  by  some  appropriate  weighting 
of  the  data. 

Another  source  of  error  in  the  measurement  of  travel  times 
is  ship  motion.   The  error  in  round  trip  travel  time  due  to 
horizontal  ship  motion  is  given  approximately  by  the  expression 

€   =  TTV  cos  e/c 

111       J\ 

where  TT   =  measured  round  trip  travel  time 

m  c 

9  =  angle  from  the  vertical  between  the  transponder  and 
ship  transducer 
Vx  =  horizontal  component  of  ship  speed  away  from  the 
transponder 
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c   =  average  sound  speed 
For  ship  speeds  of  1  knot  and  measured  travel  times  of  6  to 
12  seconds  (typical  in  deep  ocean  applications)  the  round  trip 
travel  times  will  be  in  error  by  approximately  2  to  4  msec. 
This  error  can  be  reduced  by  taking  measurements  while  stopped 
or  by  reprocessing  corrected  data  once' the  approximate  ship- 
positions  and  speeds  are  known.   Other  ship  motion  such  as 
heave,  pitch  and  roll  cannot  easily  be  accounted  for  but  when 
a  large  number  of  measurements  are  used  their  effect  can  be 
considered  as  an  uncertainty  in  transducer  depth  equivalent  to 
a  timing  error  of  approximately  1  msec. 

A  summary  of  these  errors  is  present  in  Table  1.   To  achieve 
these    errors  in  measured  travel  time  of  a  few  milliseconds 
it  is  essential  to  use  an  accurately  calibrated  system  and 
make  the  measurements  while  stopped.   These  errors  can  be 
assumed  to  have  zero  mean  for  a  large  number  of  measurements. 

To  obtain  geometric  slant  range  from  corrected  travel 
time  the  velocity  of  sound  must  be  known.   Although  the  sound 
velocity  has  spatial  and  temporal  variations,  horizontal  and 
temporal  variations  can  be  assumed  negligible.   Sound  velocity 
variation  over  depth,  however,  is  commonly  on  the  order  of 
50  meters/sec  for  5  km  deep  water.   The  structure  of  this 
variation,  the  sound  velocity  profile,  can  be  determined  from 
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Table  1.   Sources  of  errors  in  measurement  of  Travel  Time 


Source  of  Error 

(f  (msec) 

^(meters) 

Recognition  Time 
Miscellaneous  Timing  Errors 
Horizontal  Motion  (vx  < 1/2  knot) 
Other  Motion 

1.0 
1.0 
1.0 
1.0 

1.5 
1.5 
1.5 
1.5 

Total  =  {ZC2)1/2 

2.0 

3.0 
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historical  data  or  direct  measurements.   Historical  data  can 

be  expected  to  be  in  error  by  2.5  m/sec  or  more,  whereas 

present  measurement  techniques  have  an  accuracy  of  about 

.25  m/sec.   For  a  precision  navigation  system  the  use  of  direct 

measurements  is  strongly  preferred.   When  the  sound  velocity 

profile  is  known  the  travel  time  of  an  acoustic  signal  between 

two  points  can  be  determined  and  an  effective  sound  speed,  V  , 

can  be  calculated. 

Ve  =  SR/TT 

When  the  travel  time  is  measured  using  a  pulse  navigation 

system  the  effective  sound  speed  is  not  known  since  it  depends 

upon  the  positions  of  the  transponders  and  ship  transducer. 

The  slant  range,  however,  can  be  determined  approximately  by 

using  a  calculated  range  R,  equal  to  the  corrected  one-way 

travel  time  TT  multiplied  by  the  arithmetic  mean  sound 

velocity  c": 

SR  £T  R  =  cTT 
c 

ZS 

where  c~  =  /       c(z)dz 

Z.  -  ZS 


1 


i 


The  depth  of  the  transducer  ZS  is  usually  known  and  errors  of 
less  than  50  meters  in  the  estimate  of  transponder  depth,  Zj_, 
have  a  negligible  effect  on  c~.   The  arithmetic  mean  sound 
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velocity  can  therefore  be  considered  a  known  constant. 

The  ratio  of  slant  range,  SR,  to  calculated  range  R  =  cTT  , 
can  be  expressed  as  a  function  K(Ve) : 

K(Ve)  =  SR/R  =  Ve/c 
For  a  given  sound  velocity  profile  and  depths  ZS  and  Z.  the 
effective  sound  speed,  Ve,  is  a  monotonically  increasing 
function  of  the  geometric  angle  between  the  transponder  and 
the  survey  point  measured  from  the  vertical  (Vass,  1966)  (see 
Fig.  3).   Since  the  calculated  range,  R,  is  also  an  increasing 
function  of  the  geometric  angle,  K(Ve)  can  be  expressed  as  a 

function  of  R  and  approximated  by  a  second  order  polynomial 

2 

K(Ve)  =  K(R)  «  aQ  +  a-j^R  +  a„R  , 

The  coefficients,  a-,  ,  are  found  by  calculating  the  range  R 
and  the  ratio  SR/R  for  various  angles  and  applying  a  second 
order  curve  fitting  subroutine.   Travel  times  are  computed 
using  standard  ray  acoustic  techniques  (Officer,  1958) .   A 
description  of  program  ARCOR  which  performs  these  computations 
is  contained  in  Appendix  A.   The  error  in  using  the  relation- 
ship SR  =  R  is  usually  less  than  1.5  meters  but  can  be  as 

great  as  3  or  4  meters  for  angles  greater  than  80°.   When  the 

2 

relationship  SR  =  R(aQ  +  a  R  +  a~R  )  is  used  the  errors  in 

slant  range  estimates  are  less  than  .1  meters.   When  the  sound 
velocity  profile  is  accurately  known,  errors  caused  by  acoustic 
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refraction  are  negligible  compared  to  timing  and  transducer 
motion  errors. 

2)  Measurement  Errors  (Doppler  system) 

The  Doppler  system  in  its  present  configuration  (Fig.  4) 
consists  of  three  bottom-moored  CW  beacons  each  separated  by 
50  Hz,  a  phase  detector  for  each  beacon  and  a  real  time  com- 
puter processor.   Each  beacon  contains  a  crystal  oscillator 
with  a  frequency  stability  of  10"^  Hz  and  has  a  nominal  source 
level  of  166  dB  (re  ±1   Pascal  (aim).   The  beacon  signal  after 
being  processed  by  its  phase  detector  is  recorded  on  analog  mag- 
netic tape  and  is  fed  to  a  digital  computer  for  real  time  pro- 
cessing.  The  phase  detector  output  consists  of  two  base-band 
outputs  one  proportional  to  the  sine  of  the  input  phase  angle, 
the  other  proportional  to  the  cosine  of  the  phase  angle.   Phase 
variations  due  to  a  finite  receiver  velocity  v(t),  can  be  written 

as  0  (t)    =  2  TTf  tv(t)/c.   If  the  input  signal  is  represented  by 
B 

s(t)  =  A(t)  exp  (  -J  [2  TTfBt  -  (2f(t)]  | 

the  phase  detector  outputs  are  given  by 

s± (t)  =  A(t)  cos  0(t) 

and  s  (t)   =  A(t)  sin  0(t)  . 

The  instantaneous  phase  of  the  beacon  signal  can  be  determined 

from  the  ratio 

s  (t)/sn  (t)  =  tan  0(t)    ■ 
2  1 
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In  practice  the  quadrant  into  which  the  phase  angle  falls  is 
determined  simply  by  examining  the  signs  of  s-,  (t)  and  s~  (t)  . 
Phase  differences  between  phase  detector  samples  are  computed 
and  accumulated  in  quanta  of   A  /4  ^  3  cm  at  13  kHz.   The 
accumulated  phase  difference  after  T  seconds  is  of  the  form 

N 

E 

n=l 


"T  "  E   [*<V  '  0(tn-l']   .»-■«■ 


where  f   is  the  sampling  frequency  of  the  phase  detector  out- 
puts  (10  Hz)  (Porter  et  al. ,  1973). 

The  distance,  DSR,  travelled  toward  or  away  from  a  beacon 
in  time  T  is  the  accumulated  phase  0     multiplied  by  the  wave- 


length A of  the  signal 


DSR  =  0T  X 
The  error  in  measurement  of  accumulated  phase  depends  upon  the 
number  of  samples  N  =  Tf   and  the  signal  detection  probabili- 
ties  of  the  receiver.   Since  the  samples  can  be  considered  to 
be  independent  measurements  the  variance  of  the  accumulated 
phase  is  the  sum  of  the  variance  of  the  individual  samples 

2      N       2        2 

0T        n=l    0n  #n 

and  ^n2=  f(Pe) 

where  Pe   is  the  probability  of  error  in  the  n   sample.   The 

probability  of  error  is  a  function  of  the  signal-to-noise 
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ratio  at  the  phase  detector  (see  Appendix  B) .   Table  2  shows 

some  representative  values  of  Pe    in  a  . 1  sec  interval  for 

various  signal-to-noise  ratios  encountered  in  practice.   The 

associated  errors  CTi,     and  (fi,    (T  =  1  hr)  are  also  given  in 

^n       ^T 

both  number  of  quadrants  and  meters.   A  quarter  wavelength  of 
3  cm  is  assumed.   Nominal  horizontal  range  is  calculated  for 
a  source  level  of  166  dB  (re:  1  wPa)  and  a  noise  level  of 
44  dB/Hz  (re:  1  uPa)  assuming  transmission  losses  are  due  only 
to  spherical  spreading  (20  log  r)  and  attenuation  (1  dB/km) . 
Sample  calculations  are  given  in  Appendix  C. 

Random  phase  fluctuations  can  be  caused  by  multipath 
interference  and  forward  scattering.   The  root  mean  square 
phase  fluctuations  from  these  sources  has  been  estimated  to 
be  about  .13  quadrants  (Porter  et  al.(  1973).   When  compared 
with  the  errors  due  to  ambient  sea  noise  shown  in  Table  2 
they  have  little  effect  (about  5-15%  increase  in  Get    )  and 
can  be  considered  negligible. 

Slant  range  difference  is  obtained  from  accumulated  phase 
change  by  multiplication  by  a  scaling  factor:  the  wavelength 
of  the  signal.   Errors  in  wavelength  X  =  c/fR  are  caused 
primarily  by  errors  in  sound  velocity  (the  instability  of  the 
beacon  frequency  is  negligible) .   When  the  arithmetic  mean 
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TABLE  2.   Doppler  Errors  vs  Signal-to-Noise  Ratio 

i 

10  log  S/N 

P 
e 

Cnt    (1  hr) 

^T 

X 

(km) 

(dB) 

(quad) 

(quad)     (meters) 

33 

.023 

.217 

41        1.2       0 

30          .031         .253 

i 

48         1.4      3.7 

27 

.042         .293 

56         1.7      5.5 

24 

.055         .341 

i 

65         1.9      7.5 

i 

21 

.072 

.393 
.456 

75 
87 

2.2      9.2 

18 

.095 

2.6 

- 
11.1 
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sound  velocity  of  an  accurately  measured  sound  velocity  pro- 
file is  used,  the  standard  error  in  wavelength  will  usually 
be  less  than  1  meter.   This  can  be  reduced  further  by  using 
the  average  sound, 

S  =  length  of  travel  path/travel  time, 
once  the  approximate  survey  positions  are  known. 

In  summary,  when  operating  above  the  pulse  system  thresh- 
old/  slant  ranges  can  be  measured  by  the  pulse  system  with  an 
accuracy  of  3  to  5  meters.   Under  the  same  conditions  slant 
range  differences  can  be  measured  with  an  accuracy  of  1.2  to 
2.5  meters  for  measuring  intervals  of  1  hour. 

3)  Selection  of  Survey  Coordinates 

The  error  in  estimating  beacon  coordinates  also  depends 
on  the  error  magnification  terms  of  the  beacon  coordinates, 

Mi  =  4a±,       i  =  1 ,  2  ,  .  .  .  ,  6 
where  the  a-  are  the  i    diagonal  elements  of  the  A  =  (X  V  X) 
matrix.   These  terms  are  dependent  upon  the  number  and  relative 
locations  of  the  survey  points.   In  general  the  error  magnifi- 
cations will  vary  inversely  with  the  degrees  of  freedom  which 
increase  with  the  number  of  survey  points  above  a  certain 
minimum.   For  example  the  slant  range  difference  equations 
require  at  least  9  survey  points  for  a  solution  to  exist,  and 
the  degrees  of  freedom  for  this  set  of  equations  is  equal  to 
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the  number  of  additional  survey  points.   Comparison  of  survey- 
techniques  should  be  based  on  the  same  number  of  degrees  of 
freedom.   For  example  a  pulse  survey  of  7  points  (1  degree  of 
freedom)  should  be  compared  to  a  Doppler  survey  of  10  points, 
and  a  pulse  survey  of  10  points  should  be  compared  to  a  pulse- 
Doppler  survey  of  10  points. 

A  measure  of  "good"  survey  locations  is  the  trace  of  the 
beacon  coordinate  covariance  sub-matrix,  that  is,  the  sum  of 
the  squares  of  the  error  magnification  factors  of  the  beacon 

coordinates : 

6  6 

TRACE  =   £    Mi   =   E    ai- 

i=l        i=l 

The  trace  is  a  function  of  the  beacon  and  survey  location 
coordinates,  i.e.  the  parameter  vector  0,  and  will  be  denoted 
by  TR ( 9) .   For  a  given  number  of  survey  locations  there  will  be 
a  set  or  sets  of  survey  locations  for  which  the  trace  is  a  mini- 
mum.  A  typical  deep  ocean  deployment  of  the  pulse-Doppler 
navigation  system  will  be  discussed  for  illustration. 

The  most  commonly  used  3-beacon  acoustic  navigation  net 
is  ideally  an  equilateral  triangle.   The  baseline  length,  i.e. 
the  length  of  the  sides  of  the  triangle/  are  typically  on  the 
order  of  the  beacon  depths:  about  5  km.   For  this  configura- 
tion of  the  net  the  covariance  matrix  was  evaluated  for  several 
sets  of  survey  locations,  each  consisting  of  10  survey  points. 
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A  few  examples  are  given  in  Figures  5A  through  5E.   The  value 
of  the  trace  of  the  beacon  covariance  matrix  for  each  arrange- 
ment is  listed  in  Table  3.   The  trace  for  the  corresponding 
survey  with  only  pulse  data  is  given  for  comparison.   A  more 
rigorous  technique  for  determining  "optimum"  survey  locations 
was  also  used.   Powell's  algorithm  to  find  a  (local)  minimum 
of  a  function  of  several  variables  was  applied  to  a  function 
G(9)  =  TR(0)  +  C ( 9) ,  where  C (9)    is  an  arbitrary  function  added 
to  the  trace  to  constrain  the  survey  locations  within  the 
maximum  range  of  the  system  (Powell,  1965) .   No  significant 
improvement  to  the  starting  value  of  G(0)  v  3.6  was  realized 
although  several  survey  geometries  were  tried. 

4)  Computer  Simulation 

A  computer  program  called  PDSRV  was  developed  to  solve 
the  survey  problem  using  the  parameter  estimation  techniques 
discussed.   It  is  a  modification  of  a  Biomedical  Computer  Pro- 
gram (BMD  07R)  (see  Dixon,  1973).   The  program  requires  pulse 
data  from  at  least  six  survey  points  and/or  Doppler  data  from 
at  least  ten  survey  points.   In  addition  to  the  data,  the 
user  must  input  the  average  sound  velocity  and  the  acoustic 
refraction  coefficients  found  by  program  ARCOR.   The  beacon 
frequencies,  the  depth (s)  of  the  transducer  and  initial 
estimates  of  the  beacon  positions  must  also  be  provided. 
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FIGURE  5A,    TRACE  =3.6 
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TABLE  3.   Comparison  of  the  Trace  of  the  Beacon  Covariance  Matrix 
for  Pulse/Doppler  and  Pulse  Surveys 


Figure  # 

Trace 
Pulse  &  Doppler 

Trace 
Pulse  Only 

5A 

3.6 

6.8 

5B 

4.0 

8.3 

5C 

4.3 

10.4 

5D 

17.0  . 

58.2 

5E 

138 

325 
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Maximum  and  minimum  values  for  the  beacon  coordinates  are 

constrained  to  provide  more  rapid  convergence.   The  constraints 

are  set  as  follows: 

X.  =  X.  +2000  m 
3  3 

Y.  =  Y.  ±2000  m 
3  3 

Z-  =  Z.  ±1000  m 

3  3 

where  (X  •  ,  Y.,  Z.)  are  the  initial  estimates  of  the  j*-"  beacon 
coordinates.   Other  parameter  bounds  may  be  specified  if  desired. 
Although  it  is  unlikely,  even  with  the  most  rudimentary  primary 
navigation  system,  that  initial  estimate  errors  will  exceed 
those  specified  above. 

The  program  consists  of  three  main  sections:  (1)  initiali- 
zation and  data  input;  (2)  minimization  of  weighted  squared 
error;  (3)  determination  of  the ■ approximate  error  covariance 
matrix  of  the  beacon  positions.   For  convenience  the  program 
is  divided  into  7  subroutines  as  illustrated  in  Figure  6. 
A  description  of  each  subroutine  is  given  in  Table  4. 

Two  additional  subroutines  were  used  to  simulate  survey  pro- 
blems.  Subroutine  READ  calculates  data  for  a  given  set  of  survey 
locations  and  beacon  coordinates.   Subroutine  NORMA  generates 
Gaussian  distributed  numbers  with  a  specified  mean  and  standard 
deviations.   These  are  used  to  simulate  data  with  errors. 
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TABLE  4.   Description  of  Program  PDSRV 


Name  of 
Subroutine 


Description 


MAIN 


DINPT 


Reads  the  following  user  inputs: 

1.  Transducer  depths 

2 .  Convergence  criteria 

3.  Bounds  for  beacon  parameters 

4.  Beacon  coordinate  estimates 

5.  Survey  position  estimates  (optional) 

6.  Bounds  of  Survef  position  (optional) 

1.  Reads  input  data  from  specified  device (s) 

2.  Calculates  slant  ranges  (SR)  from  travel 
times 

3.  Calculates  changes  in  slant  range  (DSR) 
from  accumulated  phase  change 


MINIZ 


XPRMX 


1.  Computes  the  mean  squared  error  (EE) 
from  the  weighted  sum  of  squared  errors: 
EE  =  1/DOF)  SS(§) 

2.  Checks  that  EE  is  less  than  previous 
value;  if  not  the  step  size  is  halved 

3.  Checks  convergence  criteria 

1.  Computes  the  elements  of  the  X  V   X 
matrix 

2.  Calculates  the  weighted  sum  of  square 
errors: 

n 

ss(£)  =  2"  iyL-fL)   wi 

i=l 
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TABLE  4.   Continued 


Name  of 
Subroutine 


Description 


Calculates  values  of  the  slant  range 

and  slant  range  differences  (f-)  using 

latest  estimates  of  beacon  and  ship 

positions 

Calculates  the  elements  of  the  X  matrix, 

i.e.  the  derivatives  of  the  slant  ranges 

and  slant  range  differences: 


ID 


d 


e 


A 

9 


1.  Inverts  the  X  V  X  matrix  in  a  step-wise 
fashion 

2.  Computes  the  step  sizes  A    8 


OUTPT 


Computes  and  outputs  the  following 

1.  Final. estimates 

2.  Correlation  matrix  for  Beacon  Positions 

3.  Root  mean  square  error  (standard  error, <f  ) 

4.  Error  magnification  terms  of  the 
covariance  matrix 

5.  Residuals  and  standard  errors  of  each 
measurement  (optional) 
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Typical  output  of  a  simulated  survey  problem  is  shown  in 
Tables  5a  through  5d.   Table  5a  lists  the  beacon  and  survey 
coordinates  used  by  subroutine  READ  and  the  errors  in  data 
generated  by  subroutine  NORMA.   Table  5b  lists  the  simulated 
data  which  is  input  to  program  PDSRV.   The  final  estimates, 
correlation  matrix,  error  magnification  factors  and  standard 
errors  of  the  beacon  coordinates  are  also  shown  in  Table  5b. 
Similar  information  for  the  survey  locations  is  listed  as 
illustrated  in  Table  5c.   The  residuals, 

SR(observed)  —  SR(estimated)  and 
DSR(observed)  —  DSR(estimated) 
can  be  compared  with  the  estimates  of  standard  error  of  each 
measurement.   This  output,  labeled  SURVEY'  ERRORS,  is  optional 
and  is  shown  in  Table  5c. 

After  each  iteration  of  the  minimization  process  the  mean 
squared  error  and  latest  estimates  of  the  beacon  positions 
can  be  output  as  shown  in  Table  5d.   The  number  of  the  iteration 
and  the  number  of  times  the  step  size  was  halved  is  also 
included. 

Survey  problems  with  measurement  errors  of  more  than  750 
meters  have  been  simulated.   The  minimization  process  converged 
in  every  case  and  the  beacon  coordinate  estimates  were  nearly 
normally  distributed  about  the  true  values  with  the  standard 
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Tahle  5a.     List  of  beacon  and  survey  coordinates  used  for  survey 
simulation.     Random  normal  errors  generated  by  Subroutine  NORMA. 
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Table  5b.    (1)   Simulated  Pulse  &Doppler  measurements.      (2)   Final  estimates  of 
beacon  positions.      (3)  Correlation  matrix  for  beacon  positions. 
(4)  Error  magnification  factors  and  estimated  standard  errors  for 
each  beacon  position. 
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-Table  -5c 


(1  )-Estimates  of  survey  positions-  (XjY$  with  their  associated 
error  magnification  factors  and  estimated  standard  errors. 
(2)  Residuals  of  slant  range  measurements   (1-10)  and  slant 
range  difference  measurements   (2-10). 
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Table  f>d.     Beacon' coordinate  estimates   (X2,X3,\Y3,Z1,Z2,Z3)  at  each  iteration. 
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deviation  approximately  the  same  as  that  predicted. 

Single  simulations  with  measurement  errors  of  5  to  800 
meters  were  made  to  determine  the  stability  of  the  estimates 

A 

of  standard  error,  CT ,  and  the  approximate  error  covariance 
matrix.   Table  6  lists  some  examples  showing  the  simulated 
measurement  error,  CT  .  the  estimated  standard  error,  U    ,  and  the 
estimate  of  the  TRACE  of  the  beacon  error  covariance  matrix. 
The  stability  of  the  trace  indicates  the  linearity  of  the  model 
in  the  region  around  9^. 

The  stability  of  the  trace  has  another  important  implica- 
tion: suppose  some  "optimum"  set  of  survey  points  have  been 
determined  for  a  particular  navigation  net;  the  actual  survey 
points  used  might  differ  from  the  optimum  by  hundreds  of  meters 
due  to  navigation  errors  while  conducting  the  survey;  the 
resulting  errors  in  the  estimates  of  beacon  coordinates,  how- 
ever, would  not  be  very  different.   The  primary  geometric  factor 
which  affects  the  errors  of  the  estimates  is  the  distribution 
of  the  survey  points  not  their  exact  location. 
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TABLE    6, 


cr 

1        A 

cr 

Trace 

5 

4,8 

3.96 

50 

47,8 

3.95 

500 

513 

4.76 

600 

640 

3.92 

700 

751 

3.98 

800 

820 

4.35 
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IV  DISCUSSION 

Using  the  parameter  estimation  technique  and  the  procedures 
described,  the  relative  positions  of  acoustic  navigation  beacons 
can  be  determined  with  an  accuracy  of  ±1  meter. 

A  non-linear  model  is  developed  for  a  combined  pulse-Doppler 
navigation  system  whose  measurements  of  travel  time  and  accumula- 
ted  phase  are  perturbed  by  additive  zero-mean  Gaussian  noise. 
The  maximum  likelihood  estimation  for  the  beacon  coordinate 
is  found.   The  numerical  technique  used  to  solve  the  non-linear 
estimation  problem  is  a  modified  Gauss-Newton  method.   This  method 
consists  of  approximating  the  non-linear  model  with  a  linear 
expansion  and  iteratively  solving  the  linear  model  using  standard 
least  squares  estimation  techniques.   An  approximate  error  co- 
variance  matrix  is  found  from  which  errors  in  the  beacon  coordi- 
nate can  be  estimated.   Computer  simulation  has  shown  this 
technique  to  be  stable  for  measurement  errors  within  (and  beyond) 
the  range  of  errors  encountered  in  practice. 

Many  sources  of  error  are  present  in  the  survey  problem. 
The  preceding  sections  detail  the  several  steps  to  be  taken  in 
order  to  achieve  an  accuracy  of  ±1  meter  in  the  beacon  estimates. 
These  can  be  generalized  as  follows: 

1.  Sound  velocity  should  be  determined  with  a  precision 
instrument,  for  example,  a  Conductivity-Temperature-Depth  (CTD) 
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probe  (see  Gregg  and  Cox,  1971).   Errors  in  the  measurement  of 
sound  velocity  will  usually  bias  the  estimates  but  will  not 
affect  their  variance. 

2.  Corrections  should  be  applied  to  the  measurements  to 
account  for  acoustic  refraction.  When  these  corrections  are 
applied,  the  errors  due  to  acoustic  refraction  can  be  consi- 
dered negligible. 

3.  Travel  time  measurements  should  be  made  while  the 
surveying  platform  is  dead-in-the-water,  or  nearly  so,  to 
eliminate  errors  due  to  horizontal  platform  motion.   This  is 

a  potential  source  of  large  errors  in  the  pulse  system,  although 
it  does  not  affect  the  accuracy  of  the  Doppler  system  measurement, 

4.  The  time  interval  between  selected  Doppler  survey  points 
should  be  less  than  1  hour  since  errors  in  the  measurement  of 
accumulated  phase  increase  with  time.   Intervals  of  approximately 
1/2  hour  each  are  recommended. 

5.  Survey  locations  which  will  minimize  the  error  magnifi- 
cation factors  of  the  beacon  coordinate  estimates  should  be 
used.   The  choice  of  survey  locations  is  critical  to  obtaining 
accurate  estimates  and  can  also  affect  the  ability  of  the 
estimation  technique  to  converge  to  a  global  rather  than  a 
local  minimum.   Several  examples  of  good  survey  point  selections 
are  given  in  Section  III.  3  above. 
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6.  Pulse  and  Doppler  measurements  should  be  made  simul- 
taneously at  no  less  than  10  survey  locations.   The  use  of 
more  than  10  survey  locations  will,  in  general,  reduce  the 
errors  in  the  beacon  estimates.   For  example,  the  use  of  40 
survey  locations  can  reduce  these  errors  by  a  factor  of  1/2. 
In  most  cases  the  capacity  of  the  processing  computer  will  limit 
the  maximum  number  of  survey  locations  which  can  be  used 
practically. 
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APPENDIX  A 

DESCRIPTION  OF  PROGPAM  ARCOR 

Pulse  navigation  systems  calculate  the  distance  between 
two  points  from  the  measured  time  of  travel  of  an  acoustic 
pulse  between  those  two  points.   This  distance  can  be  approxi- 
mated by  multiplying  the  travel  time  by  an  average  sound 
velocity.   This  is  an  approximation  since  the  sound  travels 
in  a  refracted  path  rather  than  a  straight  line.   The  purpose 
of  Program  ARCOR  is  to  improve  the  approximation  by  accounting 
for  the  refraction  of  sound  in  water.   For  a  given  sound 
velocity  profile  and  the  depths  of  a  sound  source  and  receiver, 
this  program  uses  standard  ray-tracing  techniques  (Officer, 
1958)  to  find  the  time  of  travel,  TT,  of  an  acoustic  pulse 
between  the  two  points  is  the  travel  time  TT  multiplied  by  the 
arithmetic  mean  sound  velocity  c: 

R  =  CTT   t 

The  actual  distance,  SR  is  given  by  the  expression 

7  2  1/2 

SR  =  (XZ    +    Z")±/Z 

where  X  is  the  horizontal  distance  traveled  and  Z  is  the  verti- 
cal distance  between  source  and  receiver.   The  ratio  between 
the  actual  distance  and  the  approximate  distance,  K^  =  SR^/R^, 
is  calculated  for  several  travel  times.   The  set  of  approximate 
distances  and  their  ratios  (Ri,Ki),  i  =  1,N,  are  fit  to  a 
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Figure  7.   Notation  used  for  Program  ARCOR 
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second  degree  polynomial  using  standard'  least  squares  re- 
gression: 

SR 

K  =  =  a   +  a-^R  +  a2R 

R 

where  the  a.  are  the  acoustic  refraction  coefficients.   The 
actual  distance  SR  can  therefore  be  calculated  when  the 

approximate  distance  R  =  cTT  is  known:' 

2 

SR  =  R(a„  +  a,R  +  a0R  ) 

0     1      2     • 

The  set  of  approximate  distances  and  their  associated 
ratios  are  generated  by  calculating  the  necessary  parameters 
for  various  angles  9~  at  which  sound  leaves  the  source  (see 
Fig.  7) .   For  each  angle  90  there  is  a  constant  given  by 
Snell ' s  law 


sin0Q                  sinSi 
=  p  (constant)  =  


C0  Ci 

where  c~  =  sound  velocity  at  the  source  and  c.  =  sound  velo- 

city  at  interface  i.   This  constant  and  the  sound  velocity 

gradients 


gi  = 


c.  -  c.  , 
i     l-l 


zi  "  zi-l 


determine  the  travel  time  and  the  distance  an  acoustic  signal 
must  travel  to  reach  the  receiver  depth,  assuming  constant 
gradient  in  each  layer.   The  angle  at  which  sound  leaves  a 
layer  is 
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9i  =  sin    (pc)  . 
The  travel  time  for  each  layer  is 

1       tan  9-/2 

t  =  In  ( ). 

g       tan  ©i_1/2 

The  horizontal  distance  traveled  is 

1 

x.  =  (cos9.    -  cosG  )  (Officer,  1958). 

l  l  — J.        -; 

pg  1 

For  each  starting  angle  9Q  the  total  travel  time  and  horizontal 

distance  traveled  are  found  by  summing  the  above  quantities 

for  each  layer.   We  then  have 

Total  horizontal  distance:  X  =  5"*  x . 

£-i     i 

Vertical  distance:         Z  =  source  depth-receiver  depth 

Total  travel  time:        TT  =  2J  t 

2     2  1/2 
Distance  between  source  and  receiver:  SR  =   (X   +  Z  ) 

Approximate  distance:      R  =  cTT 

Ratio:  K  =  SR/R 

The  program  computes  the  first  starting  angle  by  defining  an 
angle,  c<  =   90°  and  setting  the  departure  angles  equal  to 
90°  -  o4  .   Each  successive  °^~   is  reduced  by  a  factor,  for 
example  .95,  and  a  new  departure  angle  is  computed.   The  values 
of  K  and  R  are  fitted  with  a  second  degree  polynomial  to  pro- 
duce the  coefficients,  a..   A  typical  output  is  shown  on  page 
73.    A  listing  of  the  program  begins  on  page  74.     The  pro- 
gram was  originated  by  W.M.  Marquet  of  the  Woods  Hole  Oceanographic 
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Institution  and  includes  calculations  of  coefficients  for  the 
case  of  a  receiver  at  the  same  depth  of  the  source.   This 
particular  computation  is  not  used  in  the  survey  problem  and 
is  not  discussed  here.   The  original  version  of  this  program 
was  called  Program  SETUP  and  written  by  Mary  Hunt  and  Roger 
Goldsmith  of  the  Woods  Hole  Oceanographic  Institution. 
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ARITHMETIC    MEAN    SOUND    VELOCITY    = 


1512*73^ 


?« 


.00 
7. "3 


SR 

4615 
'16  5  3 


5^ 
J  J  ^ 

i3t 

1,6^ 


5.  ^8 
5.  VI 
0.83 
1 .  53 
9.  ^0 
5.^4 


4  0  79 
47i9 
4773 
_4  R  4  0 
4919 
5009 
51  1  0 


0  0 
88 


3.U709 
3.  "7  68 


THETA 

.00 
3.60 


28 

62_ 

66 

''3_. 

16 

25 


3.  ''19  3  5 
3.1202 
3.155  9 

_3_s_£8ILL 
3.2  521 
3.3H7 


99901 
99JL21 


_R 

4  645.  (|  3 

4  65'!  .32 


99QQ1 

235L9JL 

99901 

29_?  0.1 

9  9901 

25£JJL 


4679./2 
472^.05 

4  77'»  .11' 

4  919.59 
_50_fi5^6P_ 


DEEP    TO    SHALL  O'.v  (  4650.1'^    TO 

JVSK_(__1  )_  = .9^986l6E+  0.0 

ASK(     2).    =     ".  5cl6?5  90E-i;^ 

ASK(     3)     r        .2213737E-1J 


5*00  ) 


Output  of  Program  ARCOR.  X=Horizontal  distance  between  source  and 
-receiver;"  SR  =  Slant  range ;—TT- -"Travel- time; ^  TiKTA  =  Angle  of 
departure  from  source;   K  =  SR/R;  R  =  cTT  =  approximate  slant  range. 
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XDOS    -    FORTRAN     IV    COMPILER 


00*>1 
00^2 


FTN»L 


PR 
DO 
HP 


00?  5 
00»6 

00*7 
J  0___.8 

gP»3    f; 


00^  1 

J_0{2 
0  0i  3 

JjgML 
0  0  J  5 
0  0  ■  j  6 
0  0  _  7 

"a  0 1  g 
egg 

00?  1 

jy>£2_ 

0  0  __  3 
00.^ 

0  0^  5 

'0  0£'7 

JL*!£JL 
0  &_  9 

_0_P' 
0  0  ?  i 
B0-2 
0  0?  3 
Dp  ">H 
0  0?  5 
0  0-5  6 

"06?"7 

JS  0?  8 
0  0-9 
0B'' 
H  0'J  I 
0  0_<  2 

~00'<  3 
0_0_.  'I 
00_«5 
0  >'"  6 
i •["<  7 
l'0'<O 
0  0  '*  9 


FI 
DO 
DO 
DO 
DI 

_D_I 
Dl 
DI 
DI 
DI 
EG 

JU- 
KI 
NO 


OGRA 

UBL  E 

DOS 

rjDS 

UBL  E 
UBL  E 
UBLE 
MENS 
MENS. 
MENS 
MENS 
MENS 
MENS 
U I V  A 

~SE 

-    1 

._.. ._ 


M._*R 
PRE 
VER 


COR 

c  is 

SIO 


ION    VERSION    FOR    HP 

M    BY    MARY     HUNT     &    'jAMFS     DURHAM 


AUG    7H 


RAY 
.  PRE 
PRE 

pre 

ION 
ION 
ION 
ION 
ION 
ION 
LENC 

•JLli 
T  UN 


[■  E  N 
CIS 
CIS 
CIS 
R  I  5 
DTD 
Z.»5 
GJ_5 
COE 
SAV 
I  ( 
55_2 
IT 


DING  COEFFICIENTS 
ION  PI  »X  I  »  XF  iXD  »_XM 
ION  TH02.f|-i2»TT»AR 
ION  XS*TS_CS_,__.L_2__R 
0'  »DX(5PJ 
C  <  3)  'OTSC (3) 


AX_»XC 

G'  AL  P 

T  'ASN 


)SH  . 
hi » AL 
♦  XT 


py  .t 
PhT  . 

R  A  C  T 


MGi  TH 
R  F  A  C 
»Pi  >\L 


0  )  >  C  l  5  0  ' 

R  »  »    CB(^0) 

F  «  3  i  ' ) 

E  '  3  ) 

COEF. 1,1)  *DTDC), (C 
_j_L3__  C  0 

DEVICE  CODES 


OFF 


i2!   .L  TS(  I  »  (COLF  (1  »3)  'S/  7£ 


C 

C     < 

c 


4     PAUSE    FOR    SENSK     SWITCHES 


_1 l< 
10 
1*' 


fc'RlTE  l  NO  » 1002) 
02    FORMAT  ("     PRCGf-AM     ARCOl'H     HP     VERSION"/"     SET     SFN 
1"     SSW     I     Or.'    *    S'VfcL     PROFILE    LISTED     Or;     'UTF 
3"    SSW     3     ON    "     DEEP     TO    DEEP     CALCULATION    SKifPtL 
1"     SSW     7    ON    ~     RAY     BENDK.G    CALCULATION    SKIPPED' 
L!l_SS'*  JU    ON    -    S.VEL    PKOFIL  E  .  AND.  CC-fF  FS     O'lTH     , 
WPI  TE.  NO"'  if''") 
:.i'3    FORMAT*/"     HL     INPUTS    APE    FREE    FIELD"/"     FEADY'* 
READ  <Nl'»  100'*')     IANS 

/.'(     FORMAT(M)  _        _       

WRITE.  NO  »  l?l'l  > 

1  __  FORMA  T  (_____  L  1ST    DEVICE?       -"' 

REAdTnI,*)  LP 


DO  2  I  r  li3 
COEF  U  »I  »  =  1__ 
COEF (2, I >     z     0.» 
COEF  (3,J  '  r  0." 


'/ 


U. 


ito  ■  ; 
'  / 


At':  P  "! 


/  ) 


I  i 


2  CONTINUE 

t    INPUT  SOUND  VELOCITY  PROFILE 

WRITE 'NO '10 ^5) 
_  5  F_0  R  M  A  T  _(___'_  _S  •_  V  E  L  _♦__  INPUT  LEVI  _________  .- !'  _. , 

READ  (f.I  , '  )  1  SVC 

WRITE.1  NO' 101'6) 
'06    FORMAT  ("     INPUT    PROFILE'""'     DI  f'TH  ( }    )  ,  S  .  Vi  L  !  r  /r  ' 


"/ 
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XDOS  -  FORTRAN  IV  COMPI L  ER 


0  0?  f. 
00?  1 

0  0^2 
00^3 


1"    ENTER    NEGATIVE    NUMEiER    WHEN    DONE!"' 

I_=_ 

5    CONTINUE 

1    =    1     +     1 


0  0'J'1 
J_0-_5, 

0  0'-^  6 
J_0__7_ 

0  0?  8 

0  fp  9 


READdSVC  i*'     Z  I  I  )  'C  (I) 

I_F  i  Z  (  I  )  )__  2 18  __8 

8    CONTINUE     . 


NO    MORE    THAN    5^     ALLOWED 


00f--f;  IF  5  I    "    50)    b ,  3 1: , 3 i" 

00t-l 12    CONTINUE 

0  0b 2  1=1-1 

J____C         _     .       . 

0  0f<4       C 

__t_       C     _* INITIALIZE    V  ART /PIES 

00fc6       C  SET     X     INITIAL'     TERMINAL     AND    INCREMENTAL    VALUES 

_j_j____jC 

0  0(0  3G    CONTINUE 

0JJC_9 m5v    -     I       

00/"  IF  f'iSSWli"))     31»32,22 

E|3a  31    WRITE.1  LP*  30^0) 


00.2  3^30    FORMAT  (/l"       I  Z  C"/) 

0  0/3  a  R  I TE  |_L  P_'_3  £  _  )_    (  I  »Z(  I  )  __J  I'  »I~1>NSV' 

0  0 /  <»      ~  3 * ''■:■ 1     FOR '■•  AT  (  I  u  » 2F  10  •  4  ) 
J  0£5  32     IF  M  55W  (  7  )  )     2B_I  »  3  ■';  ,  '  f| 

00/6""  34    CO  NT  I;  UE 

a  0/7 XI :    2_*_D0 

00/8  XF       ="'2000    •    D- 

___j_/9    _____  XE       =    2BJ'J5.i_D0_ 

00«..:                         WORD    =    2 
____0___1  M0RD1     =    MQRD    +    1 

S015  2  50  'CONTINUE 

B0H3  ITER    =     -1  


00t='»  C           INPUT  TYPICAL  DEEP  DEPTH 

0  0*5  WRITE  (NO____I_0^0)  

0gB6  4  M0  FORMAT  (  "  BEACON  C~EPTH?  '-"") 

£0?1Z    _     READIIIIi*)  Z? 

00^8  C~         "INPUT  SHALLOW  RECEIVER  DEPTH 

0  0t<9  WRITE l NO' 4fll ?)   


0P1'.-' 
00V1 
00V  2" 
00L'3 
00V  4 
00^5 


4^10  FORMAT t"  DEPTH  OF  RECEIVER  1?  -" > 
READ(NI»«)   ZR 


_i!lL*_SI*flJ_pE.EP  TO  DEEP  PROCESDSING 
FIND  AVERAGE  VELOCITY  AT  TYPICAL  DEEP  DTPTH. 


HP1' 6 
00^8 


71'     CONTIT  UE 

IF  l  ^Sw  ( 

75  CONTINUE 
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ARC  Op 

XPOS  -  FORTRAN  IV  COMPILER 

I 

00^9 

vim 


CPAR  =  Yp(Z0»Z»C  »NSV' 

FIND  THE  D T  EP  DEPTH  STARTING  L A_Y E R_ 

DO    110     II    =    l.NSV 

IF     (ZH    -    Z(  II  )  ;     3  3  5*3  15  »  .110 


01"  3 
01"U 
01»'5 

01ii"7 
01»8 


110  CONTINUE 

II  =  NSV 

?15  CONTINUE 


FIND  GRADIENT  OF  SOUND  VELOCITY  AT  DEEP  DEPTH 
G  (  II  )  r  l C  (  1 1  )  -  C  (  I  I  - 1  >  I  /  '  Z ( I  I)  -  Z( I  I -1  '  ) 


01  i1  9 

01  j  3 

014  l 

0.3  4  2 
0  li  3 


IF_GRA_DIENT     I  S_NEGAT_I  VE  i_  D_E  E  P_W  AT  E  R_A  N  Al  YSI  S_  1  S_  SK  IPPFD   

IF     (G'lIM     125*130.135 

_125_WRITE     lNO«9.3_0.2j .     

9J02  FORMAT! 5'H  NEGATIVE  GRADIENT  -  DEEP  TO  DEEP  PROCESSING  BYPASSED) 

GO  TO  10? 0 „.____. 


01  A  5 
03?6 
Sii  7 
Blj  8 
01J  9 

0 1-_- 1 

01^3 
0_llf4 
31*5 
Li  3  «?  6 


130    WRITE     (N0»9l  01 ' 
9 103     F0RMATU2H    GR/-CIENT_£     0    FOR  .PEEP.  TO  „DEEP  _PPCC  ESDSI  NG_/ 
+  "  20H    K     =     3  ",'c  ""FOR     ALL     X     ' 

GO    TO    10  00 


CALCULATE    [0    LOOP    PARAMETERS 


I'll  7 
01*8 
01*9 
03?:-. 

"aiTi" 

01-">2 


C 

_c 

c 

135  CONTINUE 

IE  "=  lXF  -  XI  )  /'XD  +  0".0.  IP" 
C  SET  DO  LOOP  TO  VARY  X  FOR  DEEP  TO  DEEP 

WRITE  (I  P  »QJ  r?  >   ' 
_ 9  n  •:•  5  F 0 RJ--  A  T  (  /  /  /  /  '-■  X  j  ' 'X">  IP  Xj^ t  'HAT  JH  "Z"  Lil '_ 

"do  lies'  1"  ="if  ie 

FI   £    1  

~'xc       =  xi   +"fi*xd 

IF      (XC)     3  75'  1' "■•'•'  ,375  

C  SET    THETAO    FOR    X~  =    0 

150    CONTII.UE 


03?3  DK     :    -1.0D0 

J_lf4  R(I)     =    0.0 

01? 5  GO    TO    350 

03?6       C  CALCULATE    THET_A0_ 

01?  7  "  J75    CONTINUE 
01^8  TH0    :    2.^=»CBAR/  (XC+G  <  1 1  )  ) 


01r->9  TlrJ     r    DATAN'TH*) 

J0  1'<:';  C  CALCULATE    ANGLE    TH 

01!»1  ■    *~00    CONTINUE 
Dl'<  2  TH    =    PI     "    TH0 

01?  3 C  FIND    HALF    ANGLES 

0  3,44  TH02    =    TH0/^.0D0 


i<3''5 
0  1"  6 

.3"  7 


TH*     =    TlV2.i'D& 

CALCULATE    TIME    T 
TT     z     '  DSn.tTIi? '/rCOS1  TH-)  )     /     JDSIN(TH 


/!  COS  (TH     ?»  > 
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XDOS  -  FORTRAN  IV  COMPILER 


01?  8 

J!ll?_ 

0  1?  n 

oiDi 


TT    =    DL06(TT)/G( II) 

___    CALC UL.  ATE_  APPARENT    DISTANCE    TRA  VEL  LJ T 

RACT    =    CBAR'TT 

R  (  I  )     =    R  A  C  T 


01?  2 
01^3 


01?  4 
01?5 


•550 


CALCULATE    RATIO    K 

DK    =    XC/RACT 

CONTINUE 

DX'I )     =    DK    -    1-0D0 


01?  6 

01 -'7 


00 


WRITE     (LP*9l06>     XCiDK»DX(I' 

-COLLLli'JJ  L"__ 


0  3?  8 
01?9 

0  lt> ;-; 
01b  1 
Bit-  2 
Bit  3 


FIT    POLYNOMIAL    WITH    X     AND    K     TO    GET    COEFFICIENTS 
jC_ALJL_ttCR£T  IRj  DX  lQTDCJJ  E'MORP) 

OUTPUT     CURVE    FIT     COEFFICIENTS 

DTDC  ( 1  )     '    DTDC  i  1  )     +_1_.  0      

WRITE     (LPt9l03»     Zk"  ♦  (  I  >DTDC  «  I  )  *  I     =     l'MORDl' 


Bit' 4 
0J>  6~ 

>  1 »-  8 
b  1 1  •  9 


**     PEEP    TO    SHALLOW    CALCULATIONS 


1000 


u  xi  - 
Bill 
.'1/2 
Uli  3 
wlTiJ 
l:  I  /  5 
Bl_'6 
Bl/7 
0  1  /  8~ 
01/9 
0  1*  . 
Blbl 


i  ■»  r  a 

i  152 


0  1  8  2 

0  1>  4 
01f-5 
0  3>6 
0 1 M  7 


CONTli  UE 

cs    =    r-.f-r? 

FINE     STARTING    DEPTH    AND    LAYER 

FIND    LA YEF     GRADIENT     AND     AVERAGE     SOUND    VELOCITY 

DO~  110  0    j:    IthSV 

IF     (ZK     -    Z(J.).)_11.50.»ii8_0_Lli_?0  

CONTINUE 

J    =    NSV     1     1 __ .        

COLT  I;   UL 

_K     =__J  

ZJl     =    Z« J-J ' 

CPAR    -    Yp(Z?_»Z_«C_»NSy_» 

ZJS    r    zfj) 

CUS  .=...C _<  J) . 

JS       =    J 

Z(  J)     -    1  0 

C(J)     =    CBAR 

G(  J)     -     J_C  B  AR     -    C  ( J-l'  )/  !Z(  J)     ~    ZJl) __ 

CB.'J)     r    I'CPAR     *    C.'J-1))/2.P 

GO    TO    12-0       


I  1PP 


FIND    LAYER    GRADIENT     AND     AVERAGE     VfLOCITY 
JS     =     J 


0188 

J)lb9 

_0iyi 

0 19  2 

01^3 


12  00 
1205 


1^0  8 
121  fl 


ZJl     :    Z(J-1» 

C  0 N T  1 1 . UE 

K    =    J 

IF     (Z'J)     -    2R)     1300»1300«1208 
IF     [ZJl     "    ZR)     123.0  ♦  1<-J3tf  »12-0 
CONTir.UE 


0  1  9  t< 
011'5 
01l;6 


CBAR     =     YP (ZR»Z »C »NSV > 

G(  J)     -     (C  (  J'     -    CL-AR)  /  (Z  l  J)     -    ZR) 

CPl  J)     =     "CPAR     *     C ' J)  l /?• 0       . 
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XDOS    -    FORTRAN     IV    COMPIIFR 


«1^7  K     =     1 

01?8  ZJJ     r__ZR 

B1V9  GO    TO    1250 

0  2'^  '■'•  123B    G  t  J )     -     (  C  (  J  '     -    C  (  J- 1')/  <  Z  ( J  ) 


ZJ1) 


0  2-"  1  CP'J)     =     'C(J)     H     CiJ-J))/2.f' 

J_2_2 1*5  0  CONTINUE 

02l-13  CS    r    CS    +     CD(J'  +  (Z(J»     -    Z  Jl  ) 

J2^4 K    =_K~1 ___ 

02»5  J    =    K " 

02^6  IF.lJ-D      i  30(!  »1?75«12S*0 


02»7 

.275 

0  2^8 

02^-9 

130  0 

0  2?  3 
0  2.  1 

"  1>20 

02.  2 

K24  3 

C 

15  2.  4 
"02.5 

—  ■    

0  21  6 

C 

02.  7 

0  ?  l  8 

0  2?  9 

02^  • 

920] 

02.1 

02^2 

02.3 

920  2 

:  2*- ') 

02*  5 

G2_6 
S2_  7 

C 

0  2_  8 
0  2.9 

1  1  '■  0 
C 

02-^ 

0  2"1! 

0  2  ?  2 

0  2-' 3 

02.3tJ 

0  2-5 

02:=  6 

0    DEG 


62?  7 

0258 
0  2r  9 

__2i_. 
02'J  1 

__2i'  2 
0  2,J3 
02'«  H 
:>?'•  5 


1 <<  2  0 
_«25 


l'*2  8 
1  '1 3  0 


ZJl     =    -0  •0-  t;    1 

GO_  JO    1_2J5_ 

CS    =    CS/JZ0    -ZR) 

WRITE*  LP '91j*9>_  C_S 

CONTIi  UE 

ITER    r    ITER    +    4 

SET    START  ING    ANGLE 
ALPH    -     B-.St  !' 
At  PUT    _'  PI/2.OD0 

SET  ARRAY  INDEX      

I  =  . 

RFAC  f_  •_?.5D.i 

WRITE  < NO '92*1) 

FORMATC'  ENTFR  PFCP;  N't: NT  FACTOR  BETwFFN 

READ  (NI »<  J  RFAC 

fcRlTE  !  f.O  '  Q2^!2) 

FORMAT  ("  e"nt\":'P  f-'/X  PORIZONTAt  RANGE 

x"a;._  =  3LL>..p'i-_jL:.Cj  . 

READ  IN  I »  4  )  XMAX 
WR  1TF i LP' 91  •  7) 

•   L  G I N  RAY  TRACE 
CGNTI.UE  _  

I  r^  IT  I  AL  IZF  SUPS 

X  S  _=  i '_.  0  D  2 

TS  ~    i4  .  0  [  0 

s  =  0.0r(" 

PrDSlN  IALPH'  /C  «  JS  > 
THFTA  =._ALPH*.lP.0..0p0/pI 

J  =  JS 

K  ~  J 

ZJl  _  Z (J-l  ' 

continue         _■ : _. 

CJ1  _  c(j-i' 

IF)Z(J)  "  ZR.J_..1.500t.l500.»l_/l28  

IF  i  ZJl  -  ZRJ  l'1  30'  14  "0  '  l'*'li* 

.CQNJLU  UE 

K  =  1 
ZJl  =  ?R 
CJ«  =  IP  AH 


?.     A  N 


1  -" 
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XDOS  -  FORTRAN  IV  COMPILER 


B2'J7 

0  2.'»  8 

»2M  9 


1'I4  0 


0  2VJ 
JL2£LL 

.4  2^  2 

J121L3_ 

02?  4 

02^5 


CONTINUE 

„f.IND_A  NGL„F  _B  A_Y__L .E A  VES    LAYER 
ASN    =    CJ3*P 
.J  F    (AS  N    -    1  »0DH)     Liii^jJJi^£-lIMjB 


0  2;?6 
02^7 
02^8 
02?  9 

02b 

.J2'-l 


145  0 


_LftM_ 


CONTINUE 

TH    =    A  SIN  (.  ASN  '. 

ARG    =    ASN/DSORK  l'CD*'-    -    ASN+ASN) 

T H    =    DATAN  (  , ■  r  G /_ ' 

IF  lG(  J)  )        1460»1480  .146^ 
CONTINUE 


1  470 


0  2^2 

_02j>3__ 

"i)2(-U 
02b  5 
0  2«6 

i'2»  7 


FIND     TH'E 
J1H2.    r  _ IH /  2  ,0p0 

AL?    =    ALPH/2.0D0 
..IF .'. L-l '.._]  '1 7 ?_t_lJL8 0_»  147 0 

CONTINUE 

TT    =     CTAf;(TH2XDTAN(/L2' 


TT  :     fDSlN(TH2»/DC0S.'Thf'))     /     <  DSIN  (  AL2  J  /L-COS  I AL2 

TT  =     ' -LH  O&tjT  )  )/C  (  J-'_ 

RT  t     <  "-1.0DBJ  /'PK-  (  J'  ) 

XT  t     'DCOSt  ALPH.3 -DCOS(TH)  )  *RT 

AL  r     'TH'AlPH)  *RT 

GO  T<      lur;  


)  ) 


W2&8       C  CALCULATE    LAYER    VALUES    FOR    STRAIGHT    Up     OH     C 

02^9 l^ff0_CONTINUE        • 

02/.'       C~  '  XT    ="    <Z(J"f-    ZJ4")*CTAN("lPH' 

02/3  XT    :     lZ(J)-    ZUll_iiPs_INL^Lfh)  /DCOSCALPH'  ) 

02/2"  A!      -    XT 

02'3  TT     =    LSORT     <  (Z'J)-     ZjDMZ'jl-    Z.J1  )     h     XT«^)''C! 


•  '        C. 


iLr, 


;>  2 '  '1 

~02/6~ 
0  2/7 
02/8 

i;2/9 

0  2J; 

02^2 

Ji_2£3_ 
0  2?  4 
0  2£  5 


l'«9£    CONTli.UE 

sum  x  dTsi.tance 

XS  :  XS  4  DABS1  XT) 


SUM  TR-'VEL  TIME 
TS  =  TS  *     DABSJTTJ 


SUM  ARC  LENGTH 
S_  =_S  H  AL         

RESET  ENTRY  PAY  FOP  NEXT  LAYER 
AL  PM  -  TH 

RESET  IN D  E  X  A N D  C HE C K 
K  ~     K-l 


02^6 

j  2«7 

0  2«  8 

_i)  2 !:  9_ 

B2l*f: 
02^1 
0  2V  2 
02w3 
02v'4 


3  ^95 


1  5  0  0 


J  :  K 

IF  (J  -  1)  !5  0j^»14  95^1_'i?0 

ZJl  :  -P-  0  0:-il 

GO  TO  14? 5 

"STORE  TOTALS 
CONTIf.UE 


R  A  C  T  -  X  S "  ■'  S  +  £Z0-ZRM<Z0~ZR> 
RACT  =  rSGpT(pACTi 
[<  (  1  )  -     C  S  *  T  S 
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02^5  SK     =    RACT/RlI) 

0  2^6 WR  I  TE  «  LP  '  9 1J^8L  X  5  »JR  A 

0  2^7         9108    FORMAT (E10.^»F1D.^»F 
02^8 DXJI)     -    SK     -    l.EPE 


CT»TS»THETA  »  SK  i  RJ.I  »_ 

1  E J .  't  ,  F )  0 .  2  ,  F  I  G i  .  5  »  F  I  0 


.2) 


0  2^9 
03*3 
03H1 

_0  3"  2 
0  3^3 
0  3"  4 


1*50 


0  3«5 

0  3'"1  6 

0  3*7 

JL3£8_ 

03*9 

031 


1  t  f5  0 


lb50 


i?3J_l 

ViZi  2 

"in 34  3 

_S  3J_4 

W31  5 
=  -3j  6 


0  '  i  7 
3_3L8 
0  3i  9 

13"?  I 

JL3^?L 

b  V:  3 

H  3*  ^ 

:~3i'5 
S3fT6 
0  7/7" 

■■•  .V  8 
f.«  3^  9" 

•■'.  ■»  .-< 

~03>~1 

0  3-">  2 
03>~~3 
0  334 


0  3^5 
0  3?  6 
0  3  J  7 

0  3;5  8 
03^9 

0  3'*  r 


0  3"  1 
0  3"  2 
8  3"  3 


1'/  00 

2  5 
4  i' 2  0 

I  C-^C 


2*00 

2';10 
201? 


2  '  ■  1 '( 


5 »' '-"'  0 


5>'4? 


RE 

A I  PUT 
AL  PH    = 
IF  lXS 
CONT 1 1  ■ 
Ch 
IF     (I 
I     =     I 
Fl 

coNTir; 

JLAJ.L  Jl 

OU 
DTSC  (  1 
WRITE 
IF  1 iTE 
CONTIi. 
CO  25_ 
SAVE (I 

cor.Ti- 

WRITE i 
FORMAT 

*J  COt  F_ 
READ  (N 
IF  lZR  ' 
CONTIi 
DTSC ( J 
DTSC  (2 
DTSC  (3 

coNTir: 

IF t  jSS 
CONT  If. 
IF  »Z0 

corjT  1  r. 
Z(  JS) 
C  (JS) 

contit; 

NPT  = 
WRITE  > 
FORMAT 
DO  4_5_ 
C  (  I  )  = 
WR]  TF< 
FORMAT 


INITIALIZE  ANGLE  AND  INDEX 
+_   1 ______ 

=  RFAC+ALPHT 

PI/2-0DJ?     -    AL 
-    XMAX)     155  0     » 

UE 


PHT 

1600  ♦1^00" 


DOES    NOT     flECOME    TOO    GREAT 
•'»  lf>00 


,» I  »VORD), 

COEFFICIENT 

.0 

:r  -  (  II  ,DTSC  I  I 


s 

I )  »  11  =  3  ♦  MORE  1  ) 


E(  K  THAT  INDEX 

- J18.L  14*0*140 

-'  3 

T  POLYNOMIAL  WITH  X  AND  K  COEFFICIENTS 

UE 

CRFT  (R»PX  IXS.C 

TPU1     CURVE    FIT 

)     -     DTSC.?  1. J     +. 

(L  p»9l04  '     V<  »/ 

R)    50'.l7*-fi*'\:  8 

UE 

!=•!,  1 

)     =    DTSC  '  I  )    ' 

UE  

NO' 4  020) 

(/"     DEPTH     Or     P 

!  EM  R  NE&AT 
FS  I  NOT  DESI 
I , ' )       ZP 

l95''4»3     SS»1 
Ub" 
)     =  J«8__ 

")  "="  0«'g  ' 

)_  £_p_.  r1 

UE 


ECEIV!  R    2?l'/ 

IV;      NUMEfR     IF      ANOTHER     SET     Or" 
RED'     -") 


iv  (  1  1  )  )     2i"  lr'7: 

UE 

-    Z ( JS)  )     2 Pi  4. 

UE 

r    ZJS 


2011   .2i'  14 


=    CJS 

UE 

4 

tJPT  »  5i'  ?  0  '     NSV      _       _ 

lib) 

1  =L?JiSV 

C  l  I )  / 1 0 »!  0 . 

NPT  »5^40>     Z!lHr(D 
t ; '  F  I  >'  •  4  ' 
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0  3?  4 

0  3''6 
0  3"  7 


0  348 

0  3'J  S 
0  3-1 
0  3^2 
0  3^>  3 


0  3-'  4 
J[32JL 

0  3-6 
0  3^7 
0  3^8 


'»5    CONTIflUE 

5060    FORMAT(3E20«ln? .___ 

WRITE ' NPT»5^ 60  »     (SAVE ( J>  * J  =  l «3) 

W  R ITE«NPT»506 0  »     (PTSC (J'  « J=l»3) 

WRITE  <NPT  » 5*6  0 '     (DTDC ( J » » J=l » 3 ) 

7  000    CONTINUE  

STOP 

.9j.P:3_F0RMATjy//2"H  CF.Ep  TO  DEEP    DEPTH  =  »F10»2/ 

"•         15H  ADK(*I2,3H)  =  »El<4»7)> 
g]q^  FORMAT  (///?1H  DEfP  TO  Si  ALLOW    '  ,Flfl»2»v'H  T0»FH,-2 
'  (5H  ASK  (  »I2>3H>  r'El'i.V) 

9.106    F0RMAT(3E1?.5)      

9  1 v:l    FORMA  T(///5X»1HX»C»X  f*?HSR  ,9X  »lHT  *^X  »&HT~HETA  »  5X  >  INK  »<?X 
9J.09JFORM.ATJ     /"  .ARITHMETIC  ..^EAN.SPUt  -ID    VELOCITY     =    ",     F12»4 
END 


»PH     >/ 


HiR' 


*  4       NO    ERR! 
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M_I 

MC 
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A(  I 
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A     I 
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C  (  J 
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DO 
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DO 
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+ 
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Xj_Y__ 
TT    ~ 

N    VE 
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CiN.MCURVj 

PACKARD 
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12TH    MAY     1  9  6'* 


N   .X  I  1  )  »  Y  ( 
N     A  3  5  *  fc  )  i 

FF     INDCTR 


MARCH_  l9i__ 

1)  ,C(1' 

AP  (b  t6  > ' 

•     L  =1     IS    SPCL    ROUTE     INDCTR    FOR    Y=CU>     ONLY 


=     M 
S    ORD 
IS    J 

0  ALL 
2  1  =  1 
>=g.g 
2  J  =  _ 
tjj.=0 

1  »J)  = 

*v.  - 

s  cor 

l  A    PR 

y  i  j = i 

R  =  - 


i 

ER.-Q.E_.COE 
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ARR/-YS 

5  M 
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Y    C 


lENT    MATRIX 

OL  UMN 


*    I" 

n 
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IME) 


IEMT 

IS     A 


»M     ARt--     CO 
.1 


MATR 

ux  n 

EFF  I 


IX 

IARY    MATRIX _ 

C  IENTS    OF     f.lH 


J)  »     ANC     CF    TIE    INI 


:p  .   \j 


iLE 


0  0:">5 

0  g;'  6 

00?  7" 

0  0  =  8 
0  0?  9 

0  01*  . 


19 


Y.  0_  1 
0  0"  2 
0  0"  3 
0  01*  4 
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.21 


2? 
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DO 

A(  I 

DO 

At  * 

DO 

IXP 

DO 

A(l 

THI 

DO 

AP' 

D0_ 

AP< 
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1  8 
i  J 
NO 
19 
R  = 
19 
=  J 
19 
i  J 
2? 
»M 
21 
=J 
21 


J= 

JJ 
)"=/' 

T  I 
1  = 

INC 
J- 

jJN 

J  J 
I  _A 

JJ 

C)  : 
1  = 

-1 


I  NCR 

EjL'.M 

1  I  1  J  )  ■•  X  I  J  J  )  +  1 

NCI  UDE  _J  =  1     SI 

?  » M 

R+l 

1  »M 

CR 


IXP 

fJCE 


ANY    ZERO    TO     ;>E     2...  PO    Mi  1     ALlGwEC 
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1  I  i  J)+X  I  JJ)  ** 
=  liN 

A(1*_MC  J+Y.l  JJ) 

2  » M  " 


IXP 


Jjrl »n 

,"C ) =A ( I 'MC >+Y ' JJ' 
S    FILLS     ARRAY.     NOW 
2  2    1  =  1  »M 
I ♦ 1  ) _A  (  I  »  1  ) 

2  3    Jr.-'yMC 

It  J) 


+X I JJ' ++I XP 
SOLVE    FOR    CI  J' 


[   Y     CROUT     REDUCTION 


:     A  '  1 ,  J  )  /  A  <  1  ♦  1 ) 
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00?8 

DO    28    J=^»MC 

0  0'' 9 

DO    28    T=2.M 

0  0?  0 

Sr0.0 

0tf^L_ 

TF     (I- J)     ?fi».?1',?'l 

00^2 

2M    MA    r    J-J 

00^3 

DO    2  5    K=3  »MA 

0  0"-">  '1 

25    S  =  5+AP (  I  »K  )  ♦ AP  «K iJ) 

0  0  b  5 

API  I »J)     =     A  I  I ,j)     -    S 

0  0b  6 

GO    TO    2b 

0  0b  7 

2ft    MA    r     1-1 

0  0?  8 

DO    27    Krl  »MA 

fA  0b  9 

27    S=S+AP(I»K) * AP'K »J) 

0  0b!. 

AP'I'J)     =     (A (I  'J)     -S)/AP(I»IJ 

0  0b  1 

28    CONTINUE    • 

0  0^  2 

C  (M)rAP(f--  iMC) 

gflb3 

DO    30..J  =  2»M 

0  0f-  *4 

5r0.0 

0  0'fc  5 

MArf.  -1+2 

00fc6 

MF=MA-1 

.  .jPbJl? 

DO    2  0    K=M.A  it.. 

... 

3  0C  8 

29    SrS+  AP  C-'fi.K  '  *  C.  j  K  ) 

0  0f>9 

3H    C  (MB)  =     AP(MB»f  C)     "    S 

[i0/ 

RETURN 

U  0  '  1 

fnt 

*  4        NO 

ERFO;  5  + 
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00*1 
0  0B2 
00  *>  3 

00^ 

00*5 
00»6 


YP'l 

INTERPOL  ATlOiN    FUNCTION 

WRITT 

FUNCT 

DH' EN 

YO  I  XC 


en  by  hank  perkins 
Ion  yp.1  xp'Xi y.n' 


SI  ON    X'D'Yin 

'XI  ' Y 1  'X2'Y?'      -     ( Y24  (XO-Xl) ~Y1 «  (X0-X2)  )/  ( 


XP-Xl  ) 


00V  7 

00"  8 

10 

IF«N- 

CONTI 

0  0^9 
0  01  fi 

YP    = 
RETUR 

00±1 

00]  2 

20 

CONTI 
IF'XP 

001  3 
0  01  4 

30 

CONTI 
YP    r 

0  01  5 
00  j  6 
0  0i  7 

Pi^l8 


(40 

5  0 


C0^9 

0  0*1 

0C!^2 

15*0*3 


k  £  t"  6 
5  0*7 
00*8 
0  0^9 


RETUR 
CONTI 

IF  !  XP 
COf.'Tl 
YP  t 
RETUR 
CONTI 
DO  70 
IF  <  XP 
CONTI 
CONTI 
J  =_  I 
YP  = 
RETUR 
END 


Y  0  (XP»X(H  >Y  [  lJ.'X  J[ZJ  _♦  Y.i2j. .) 

n 

t.UE  _  

-X  (N)  )        60  »  5:'  »5'-: 

:  :U  F . 

YO(XPtX(N-l)»Y(N-l'  iX'N)  ♦ Y ( N >  ) 

N  

i.UE 

I  =  1 1 N 

-  X  (  I  )  )        8 >'  »    :     i  7 

UJE 

:.UE 

Y G ( XP , X ( J" 1 )  i Y ( J- 1 '  » X  '  J  )  » Y ( J )  ) 
N 


*  +       NO     ERRO;  S* 
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APPENDIX  B 

DETECTION  PROBABILITIES  OF  WOODS  HOLE  OCEANOGPAPHIC  INSTITUTION 
CYCLE  COUNTER 

The  detection  probabilities  of  the  Woods  Hole  Oceanographic 
Institution  cycle  counter  can  be  calculated  by  considering  the 
problem  of  detecting  a  sine  wave  in  the  presence  of  Gaussian 
noise.   This  treatment  is  sufficient  because  the  estimate  of 
phase  is  not  based  on  the  optimum  estimator 

-1  S2(t) 
0   =  tan  x  

sx(t) 
where  s-i  (t)  and  S2  (t)  are  the  quadrature  components  of  the 
input  to  the  phase  detector;  the  estimate  of  phase  is  based 
on  the  signs  of  s, (t)  and  s2 (t) .   The  detection  probabilities 
of  the  sine  component  will  be  developed  for  illustration. 
The  corresponding  probabilities  of  the  cosine  component  can 
be  developed  in  a  similar  manner  with  identical  results. 
Consider  a  sine  wave  in  the  absence  of  noise 
s(t)  =  A(t)  sin  0(t) . ■ 
The  amplitude  A(t)  varies  slowly  compared  with  the  phase  0  (t) 
and  can  be  considered  a  constant.   The  phase  can  be  assumed 
uniformly  distributed 

p(0)    =    1/2  tt  ,  \0\    <  TT 
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The  probability  density  of  the  signal  p(s)  is  therefore  given 
by  the  expression: 


P(s)  = 


TTA 


1- 


A 


<  1   . 


The  phase  detector  output  in  the  presence  of  noise  will  be 

denoted  by  z(t)  =  s(t)  +  n(t)  where  the  noise  is  normally 

distributed: 

1      -n2/2  Cn2 
p(n)  -    e  (f    -   r.m.s.  noise  voltage. 

V5ifcrn 


The  probability  density  p(z)  may  be  found  by  the  convolu- 
tion of  p(s)  and  p(n)  assuming  statistical  independence  of  the 
two  processes.   This  leads  to  the  following  integral  for  p(z): 


p(z)  = 


r-r?/2C2         1 
n 


-\/2TF(/  J 

n  -' 


TTA 


1  -  ^1 


ir- 


Vo 


dn 


This  integral  has  several  representations  (se'e  Ol'shevskii, 
1967)  but  is  approximately  a  constant  near  z  =  0  for  A/fcJL >  *  L« 

Figure  8a  illustrates  the  shape  of  this  probability  density 

?  A2 

as  a  function  of  the  parameter  q^  : 


2^2  • 


The  probability  of  error  is  a  weighted  sum  of  the  proba- 
bilities of  false  alarm  P(FA)  and  false  dismissal  P(FD).   The 
cycle  counter  decides  that  the  signal  is  positive  when  the 
voltage  exceeds  some  threshold  V™  >  0..  Similarly  the  decision 
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that  the  signal  is  negative  occurs  when  the  voltage  is  less 
than  some  negative  threshold  -VT<0.   A  "false  alarm"  occurs 
when  the  signal  is  below  the  threshold  and  the  signal  plus 
noise  exceeds  the  thresholds: 

|z(t)|  >  VT  given   |s(t)|  <  VT  . 
The  probability  of  this  event,  P(FA),  is  approximately  the  same 
as  the  probability  that  the  noise  will  exceed  the  thresholds 
in  the  absence  of  signal.   This  is  represented  by  the  shaded 
area  in  Figure  8b  and  is  given  by  the  expression 

f  f     i 

P(FA)  -Z.     2  /   p(n)dn  =  2   I    — —  exp  (-x2/2)  dx 
Ju  yJ      V27T 

Vr  T4fn 

A  false  dismissal  occurs  when  the  signal  is  above  the 

thresholds  but  the  signal  plus  noise  is  below  the  thresholds: 


z(t) 


<  VT  given 


s(t) 


>  VT.   The  probability  of  this  event 
P(FD),  is  approximately  the  same  as  the  probability  that  the 
signal  plus  noise  is  below  the  thresholds  for  any  signal.   This 
is  represented  by  the  shaded  area  in  Figure  8c  and  given 
by  the  expression 

P(FD)  £  /      p(z)3z. 

-vm 

T 
In  the  region  of  interest  p(z)  can  be  approximated  by  a  constant 


p(z,s  _fA/^rn»i.  v     <     i 
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Figure  8.   Probability  Densities  of  Signal  Plus  Noise 
p(z)  and  noise  p(n). 
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The  probability  of  false  dismissal  can  then  be  found 

P(FD)  ~   VT/ 


o  2  (vT/crn) 


^a    Tr  (A/cf  ) 

n' 

It  can  be  seen  that  an  increase  in  the  threshold  voltage- 
to-noise  ratio  VT/  Cf       will  increase  the  probability  of  false 
dismissal  and  decrease  the  probability  of  false  alarm  for  a 
given  signal-to-noise  ratio.   Since  the  losses  assigned  to 
these  errors  are  approximately  equal,  the  optimum  threshold 
setting  for  a  given  signal-to-noise  ratio  will  satisfy  the 
following  equation: 

P(FA) 


=  1 


P(FD) 
or  equivalently   ^ 

I       ,    vfe  ~x2/2         vT/crn 

I    1/[2TT)  e    •  dx  =  -       (10) 

The  expression  on  the  left  hand  side  of  equation  (10)  is  a 
tabulated  function  of  Vrp/Cl  and  can  be  solved  by  trial-and- 

error  for  various  values  of  the  signal  power-to-noise  ratio 

2 

S/N  =  A2/2  CT        (see  Table  7). 

The  probability  of  error  P  in  determining  the  sign  of  the 

signal  is  given  by  the  following  expression  for  the  optimum 
threshold  setting: 

Pe  =  1/2  [P(FD)  +  P(FA)j  =  P(FD).  (11) 
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TABLE    7 


10    log    S/N 

vVcf 

/    n 

P 
e 

q;2 

12 

1.40 

.158 

.366 

18 

1.67 

.095 

.208 

24 

1.92 

.055 

.116 

30 

2.16 

.031 

.064 
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An  error  in  the  cycle  count,  i.e.  the  number  of  quadrants 
of  phase  change  between  two  successive  cycles  can  be  caused  by 
an  error  in  either  or  both  of  the  quadrature  components.   When 
an  error  (false  alarm  or  false  dismissal)  occurs  in  only  one 
component  the  cycle  count  will  be  in  error  by  ±1  quadrant. 
When  an  error  occurs  in  both  components  the  cycle  count  will  be 
in  error  by  ±2  quadrants.   The  following  probability  of  cycle 
count  error  can  be  calculated  from  the  probability  of  error  P 
previously  determined  by  equation  (11)  ". 

P(e=0)  =  (1-Pe)2 

P(e=±l)  =  2Pe(l-Pe) 

2 
P(e=±2)  =  Pe   s 

The  variance  in  the  n  "  sample  can  be  calculated  as  follows 

%  2    =  P(e=0)  (0)2  +  P(e=±l)  (l)2  +  P(e=±2)  (2)2 
^n 

=  2Pe(l-Pe)  +  4Pe2  =  2Pe(l+Pe). 
The  variance  of  error  in  the  measurement  of  accumulated  phase 
change  0"T  after  T  seconds  is 


%  2  =  Z    %  2  =  h  cr0  2      n  =  Tfs 

T    n=l    n       ^n 
where  fc  is  the  sampling  frequency.   For  a  sampling  frequency 
of  10  Hz  and  an  interval  of  1  hour  =3600  seconds  the  standard 
error  in  the  measurement  of  0T  is  (7.     z      190  CCt    . 


■91- 


Table  7  lists  the  values  of  P  ,  &h  '     and  the  optimum 

e     v>n 

threshold-to-noise  ratios  for  various  signal-to-noise  ratios 
typically  encountered  in  practice.   The  variance  is  calculated 
assuming  that  the  bias  is  zero  or  negligible.   This  assumption 
is  valid  when  the  velocity  relative  to  a  beacon  is  small. 
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APPENDIX  C 

SAMPLE  CALCULATION  OF  SIGNAL-TO-NOISE  RATIO 

In  general  the  signal-to-noise  level  will  decrease  as  the 
range  between  a  sound  source  (e.g.   CW  beacon)  and  a  receiver 
(e.g.  ship  transducer)  increases.   A  sample  calculation  is 
made  to  demonstrate  the  realtionship  between  signal-to-noise 
ratio  at  the  source  (S/N)   the  range  between  source  and  receiver 
R,  and  signal-to-noise  ratio  at  the  receiver,  (S/N)  .   The 
signal  level  of  the  CW  beacons  in  the  present  configuration  of 
the  pulse-Doppler  system  is  166  dB  (re:  l^Pa).   A  typical 
ambient  sea  noise  level  in  sea  state  3  is  44  dB/Hz  (re:  1  u   Pa) 
or  54  dB/10  Hz  (see  Urick,  1967).   The  signal-to-noise  level  in 
a  10  Hz  band  is  therefore  112  dB.   Transmission  losses  are 
assumed  due  to  spherical  spreading  and  attenuation.   The  loss 
due  to  spherical  spreading  is  20  log  R  where  R  is  the  range  in 
meters.   The  attenuation  of  a  13  kHz  signal  is  approximately 
1  dB/km.   The  signal-to-noise  level  at  the  receiver  (S/N)r  is 
therefore:  10  log(S/N)r  =  112  dB  -  20  log  R  -  R  x  10~3  .   For 
a  horizontal  range  of  10  km  and  water  depth  of  5  km  the  range 
between  a  bottom-moored  source  and  a  surface  receiver  is 
approximately 
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R  =  (52  +  102)1/2  x  103  =  11.2  x  103  m 


The  (S/N)   in  decibels  is  therefore 


10  log (S/N) r  =  112  -  81.0  -  11.2  =  19.8  dB. 
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